Median percent change: A robust alternative for assessing temporal trends

Cancer Epidemiology 37 (2013) 843–849

Contents lists available at ScienceDirect

Cancer Epidemiology The International Journal of Cancer Epidemiology, Detection, and Prevention journal homepage: www.cancerepidemiology.net

Median percent change: A robust alternative for assessing temporal trends Marco Geraci a,*, Robert D. Alston b, Jillian M. Birch b a b

University College London, London, United Kingdom Cancer Research UK Paediatric and Familial Cancer Research Group, University of Manchester, Manchester, United Kingdom

A R T I C L E I N F O

A B S T R A C T

Article history: Received 12 April 2013 Received in revised form 10 July 2013 Accepted 2 August 2013 Available online 7 September 2013

A typical summary statistic for temporal trends is the average percent change (APC). The APC is estimated by using a generalized linear model, usually under the assumption of linearity on the logarithmic scale. A serious limitation of least-squares type estimators is their sensitivity to outliers. The goal of this study is twofold: firstly, we propose a robust and easy-to-compute measure of the temporal trend based on the median of the rates (median percent change – MPC), rather than their mean, under the hypothesis of constant relative change; secondly, we investigate the performance of several models for estimating the rate of change when some of the most common model assumptions are violated. We provide some guidance on the practices of the estimation of temporal trends when using different models under different circumstances. The robustness property of the median is assessed in a simulation study, which shows that the MPC provides strong reductions in estimation bias and variance in presence of outliers. We also demonstrate how a mathematical property of the median helps addressing the issue of zero counts when estimating trends on the log-scale. Finally, we analyzed an English cancer registration dataset to illustrate the proposed method. We believe that, as a good practice, both APC and MPC should be presented when sensitivity issues arise. ß 2013 Elsevier Ltd. All rights reserved.

Keywords: Cancer Outliers Quantile regression Robust

1. Introduction In epidemiological studies, trends in incidence and mortality are of special interest. Rates are usually calculated for a prespecified number of time intervals within the study period and then plotted in order to visually assess their behavior over time. A very popular statistic used to summarize trends is the average percent change (APC). This is the percentage at which rates change between two consecutive time intervals and it is often assumed to be constant throughout the entire time period. Without loss of generality, in this paper we refer to one-year time intervals although similar arguments can be extended to the case of time intervals of different duration. Suppose we observe nt events (e.g., incident cases or deaths) in the population Pt at a given time t, t = 1, . . ., T. The (crude) incidence or mortality rate will be given by rt = nt/ Pt. The assumption of constant relative change is suitable for using a linear regression model (LM) of the type

estimated, and et is a normally distributed error term with zero mean and constant variance. Additional features of the data, such as heteroscedasticity and residual correlation, can be accounted for by using appropriate methods. Model (1) is based on the assumption that the rates vary continuously over the positive real line. This assumption has been widely used in several modeling approaches (see for example [1,2]). In effect, incidence and mortality rates, as well as survival rates, can be considered as a measure of continuous underlying risk processes. However, in order to ensure that the logarithm is defined, usually a small constant is added to observations with zero events (zero counts). Such adjustment may be source of bias, especially when the mean counts are small [3]. The regression coefficients a and b in model (1) can be estimated by solving the least squares (LS) problem min

2

T X 2 ðlog r t  a  btÞ :

(2)

ða;bÞ 2 R t¼1

log r t ¼ a þ bt þ et ;

(1)

where log rt denotes the natural logarithm of the rate observed in year t, a and b are unknown fixed regression coefficients to be

The derivative of the expected value of log rt given t, with respect to t, provides the rate at which, on average, the outcome varies for a unit change of t (e.g., one year), that is

@ Meanðlog r t jtÞ ¼ b: @t * Corresponding author at: Centre for Paediatric Epidemiology and Biostatistics, MRC Centre of Epidemiology for Child Health, UCL Institute of Child Health, 30 Guilford Street, London WC1N 1EH, United Kingdom. Tel.: þ44 2079052341. E-mail address: [email protected] (M. Geraci). 1877-7821/$ – see front matter ß 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.canep.2013.08.002

(3)

An estimate of the APC is obtained simply with APC ¼ 100  fexpðbˆLS Þ  1g;

(4)

844

M. Geraci et al. / Cancer Epidemiology 37 (2013) 843–849

where bˆLS is the least square estimate of b. Unless the distribution of the rates is truly log-normal, the APC calculated in (4) represents an approximation of the actual average relative change since, in general, Mean(log rt|t) 6¼ log Mean(rt|t). Such approximation may or may not be satisfactory in some situations. An alternative way of calculating the APC is to model the counts instead of the log-rates and to estimate the slope b by using a Poisson or log-linear regression model (GLM) of the type log mt = log Pt + a + bt, where mt is the expected number of events at time t and log Pt is an offset. The advantage of this approach is that, if the count at time t truly follows a Poisson process with mean mt, there is no need for using a LM, whose performance depends on several approximations and, ultimately, on the distribution of the rates. The Poisson model, also, can naturally deal with zero counts. However, it is not straightforward to use a model for counts when the goal is to estimate the APC of standardized rates. Another possible, well-known limitation of Poisson models is due to its parameterization, which establishes equality of mean and variance. Deviations from the Poisson assumptions can be addressed with alternative distributions for counts. For example, overdispersion (i.e., when the variability of the outcome exceeds its average) can be accommodated with a negative binomial link in GLMs. Each model has its advantages and disadvantages. In the following, we propose a robust and easy-to-compute measure of the percent change based on the median of the rates, rather than their mean, which addresses some of the limitations of GLMs, namely their sensitivity to outliers and to adjustments for zero counts in particular, and to distributional assumptions in general. Secondly, we investigate the performance of several models for assessing time trends under a variety of simulated scenarios. The models considered here include the Gaussian, the Poisson, the negative binomial (NB), and the median regression (QR). We then present an application using data from cancer registrations in England between 1990 and 2006. We conclude with some final remarks. 2. Median percent change It is well known that the LS estimator is sensitive to the presence of influential observations such as outliers with a substantial leverage. In this case, the parameter’s estimate is subject to inefficiency (i.e., high estimation variability) and bias that eventually will result in misleading conclusions. We propose a simple way to estimate the annual percent change by using a robust procedure in the sense of Huber [4] and we call it median percent change (MPC). We begin by solving the minimization problem min

2

T X jlog r t  a  btj;

(5)

ða;bÞ 2 R t¼1

which provides the least absolute deviations (LAD) estimate of b, bˆLAD . Again, we assume that the log-rates are continuous. The minimization problem in expression (5) is equivalent to finding a median regression line for the data points (t, log rt), that is Med(log rt|t) = a + bt. Median regression, in turn, can be considered as a particular case of quantile regression [5,6]. Commands for quantile regression estimation are available in popular statistical software, including R [7–9], SAS (QUANTREG procedure) and Stata (qreg command). The interpretation of the parameter b in expression (5) is analogous to that seen in expression (3). However, b is now the slope of the conditional median, that is

@ Medðlog r t jtÞ ¼ b: @t

Some of the most attractive properties of median and, more in general, quantile regression models include: (a) they do not depend on distributional assumptions; (b) they are robust to outliers; (c) they possess the property of equivariance to monotone transformations. Properties (a) and (b) are well-known. Property (c) has several important, practical implications. First and foremost, it allows us to derive mathematically an expression for the MPC. In fact, we write

Medðlog r t jtÞ ¼ log Medðr t jtÞ i.e., the median of the logarithm of the rates is equal to the logarithm of the median of the rates, or equivalently Medðr t jtÞ ¼ expfMedðlog r t jtÞg:

(6)

We noted before that this property, except in certain specific cases, does not apply to the expected value. Now, let us calculate the relative change at two time points, t and t + dt, given by

Medðr tþdt jt þ dtÞ  Medðr t jtÞ Medðr t jtÞ and recall that we assumed the model Med(log rt|t) = a + bt. By using expression (6) we obtain

expfa þ bðt þ dtÞg  expða þ btÞ expða þ btÞ

¼ expfa þ bðt þ dtÞa  btg 1 ¼ expðb  dtÞ  1:

That is, the relative change depends on the rate b and is proportional to the time interval dt. The MPC for a unitary, say one-year, time interval is simply given by MPC ¼ 100  fexpðbˆLAD Þ  1g:

(7)

where bˆLAD is the median slope obtained by solving (5). Another advantage of the equivariance property is that the median is insensitive to the constant used to replace zero counts. In mean regression, the extent of the impact of this adjustment on the parameters’ estimation will depend not only on the proportion of zeros but on the value of the constant itself [3]. There has been some discussion on the ‘optimal’ addition to make and various justifications have been given. All arguments, however, remain ad hoc and based on specific models [e.g. 10–12], with solutions whose performance needs to be assessed case by case. In contrast, median regression is invariant to censoring from below up to the median [13]. This property of quantiles has been widely exploited in regression modeling of censored data for example [14,15]. However, when the proportion of zeros is substantial, one should expect to incur considerable bias when modeling the rates rather than the counts themselves, regardless of whether one uses the mean or the median. Using a finite mixture distribution, such as a zero-inflated Poisson distribution, will in general provide a better fit. It is interesting to note that quantile estimation gives accurate results even in presence of a limited degree of zero-inflation [13]. To illustrate this point, we briefly anticipate part of the results of the English cancer registry data analysis in Section 4. Fig. 1 shows the time plot of the log-rates of oligodendroglioma (ODG) between 1990 and 2006. In 1990, there were no registered ODG cases among people aged 25–29 years. Using three different constants (0.1, 0.5, and 1), in place of the zero count in 1990, produced different estimates of the APC. This adjustment, however, did not affect the MPC. The R and Stata code to calculate the MPC, along with the data shown in Fig. 1, are provided as Inline Supplementary Computer Codes S1 and S2.

M. Geraci et al. / Cancer Epidemiology 37 (2013) 843–849

Log−rate

−3

−2

−1

The cases were generated under the processes (a–c) for each of the 9 combinations of the settings in (i) and (ii), thus giving 27 possible scenarios in total. For each scenario, 5000 datasets were replicated independently. Finally, LM, GLM and QR were fitted to each dataset in all scenarios. NB regressions were only carried out with overdispersed data. A similar analysis was conducted to assess the observed type I error rate (rejection rate – RR) at the nominal 5% level in all scenarios. This was assessed as the number of times that the true, null hypothesis of no change (b = 0) was rejected at the 5% level. For LM and QR, we introduced the customary, albeit arbitrary, approximation of adding the constant 0.5 when zero counts were generated. Each model was evaluated in terms of relative absolute error (RAE), standard deviation (SD) and mean square error (MSE) of the estimates. In addition, the average standard error (ASE), power, average length (AL) and coverage probability of 95% confidence intervals were calculated.

−4

1

−5

0.5

3.2. Standard error estimation for MPC

−6

0.1 1990

845

1995

2000

2005

Year

Fig. 1. English cancer registry data, 1990–2006. Time plot of the log-rates of oligodendroglioma in people aged 25–29 years with regression lines superimposed (linear model, dashed; median regression model, solid) which shows the effect of replacing a zero count with three different constants (filled triangles).

Inline Supplementary Computer Code S1 can be found online at http://dx.doi.org/10.1016/j.canep.2013.08.002. Inline Supplementary Computer Code S2 can be found online at http://dx.doi.org/10.1016/j.canep.2013.08.002. 3. Simulation study 3.1. Settings We conducted a Monte Carlo study to assess the impact that some common deviations from the hypothesized model can have on the inferential process. These were outlier contamination and overdispersion. In addition, we considered the influence of varying time interval lengths and different magnitudes of the expected number of events. All calculations were performed using R [7]. We generated cases using pseudo-random processes with the following settings: i. Three time intervals of 8, 15 and 30 years respectively. The annual percent change for each time period was set so as to determine a doubling of the rates from the start to the end of the period. This resulted in three APC values (steep, medium and shallow): 8.7% (8 years), 4.6% (15 years) and 2.3% (30 years); ii. Three different mean parameters (small, medium and large). At the intervals’ midpoints, these were approximately equal to 5, 170, 2060 (8 years interval); 5, 160, 1940 (15 years); 5, 150, 1870 (30 years). The case-generating processes were: a. Poisson process; b. Poisson process with outliers. These were generated by multiplying the number of cases by a factor uniformly drawn between 0.5 and 1.5. Each observation had a 10% probability of being contaminated; c. Negative binomial process with mean m and variance m + m2/n. We set the dispersion parameter to n = 2m to obtain a 50% variance inflation (heterogeneous Poisson intensity parameter);

There are several methods to calculate standard errors of quantile regression coefficients. These are implemented in R [8]: ‘iid’, which assumes independent and identically distributed errors; ‘nid’, which assumes independent but not identically distributed errors; ‘rank’, based on the inversion of a rank test either under the iid or nid assumption; ‘kernel’, based on the kernel density approach [16]; ‘bootstrap’, of which there are several variants. Simulation-based studies on the performance of most of these methods can be found in several works [6,17–20]. However, all these authors discuss coverage probabilities for sample sizes larger than those used in our simulation. We therefore conducted a more extensive (10,000 replicated datasets) prior study to choose an optimal standard error calculation method for the MPC. The type I error rate was evaluated for all methods described above, including the bootstrap method in all its variants with varying number of bootstrap replications (50 and 200). In addition, we evaluated the likelihood ratio test (LRT) based on the logistic model [19]. We focused our attention on the Poisson processes (a) and (b) as described in Section 3.1. For the uncontaminated Poisson, the ‘nid’ method performed quite well at T = 8, with RRs close to the nominal 5%. The RRs produced by the ‘iid’ and ‘rank’ methods ranged between 3 and 5%. In contrast, none of the other methods seemed to provide reasonable error rates. For an instance, the kernel and bootstrap methods were all very conservative (i.e., standard errors were overestimated leading to RRs lower than the nominal rate); on the contrary, the LRT was on the liberal side (i.e., standard errors were underestimated leading to RRs higher than the nominal rate). At larger number of years, the rank-based method showed RRs between 3.5 and 6%. These values were insensitive to outlier contamination [20]. For T = 30, the logistic LRT also provided error rates close to the nominal value. On the basis of such results, we decided to use the rank-based method for T = 15 and T = 30 in all subsequent simulations. The ‘nid’ option for T = 8 was motivated by the consideration that, being the variance of the log-rates in each year inversely proportional to the mean number of cases in that year [1], heteroscedasticity is more evident for rapidly increasing (log-) rates and, thus, for higher values of the APC. 3.3. Results For the sake of brevity, only selected results will be shown here. The complete results are reported in Inline Supplementary Tables S1 and S2.

M. Geraci et al. / Cancer Epidemiology 37 (2013) 843–849

846

Table 1 Performance statistics of the slope estimators under Poisson distributions. Abbreviations: APC, average percent change; RAE, relative absolute error; SD, standard deviation; MSE, mean squared error; ASE, average standard error; AL, confidence interval average length, RR, rejection rate; LM, linear regression; GLM, Poisson regression; QR, median regression. APC/mean

RAE

SD

MSE

ASE

Power

AL

Coverage

RR

Steep/small

LM GLM QR

5.114 4.473 5.302

0.569 0.491 0.605

0.324 0.241 0.366

0.531 0.486 0.585

0.051 0.050 0.053

2.599 2.381 2.863

0.952 0.989 0.949

0.053 0.054 0.054

Steep/large

LM GLM QR

0.228 0.228 0.265

0.025 0.025 0.029

0.001 0.001 0.001

0.023 0.024 0.025

0.860 0.954 0.770

0.111 0.116 0.125

0.944 0.981 0.944

0.049 0.052 0.056

Medium/medium

LM GLM QR

1.141 1.137 1.366

0.067 0.066 0.080

0.004 0.004 0.006

0.065 0.066 0.086

0.096 0.107 0.062

0.282 0.287 0.372

0.949 0.966 0.962

0.047 0.045 0.037

Medium/large

LM GLM QR

0.331 0.331 0.401

0.019 0.019 0.023

0.000 0.000 0.001

0.019 0.019 0.024

0.616 0.678 0.417

0.081 0.082 0.105

0.950 0.968 0.958

0.052 0.056 0.041

Shallow/large

LM GLM QR

0.483 0.483 0.609

0.014 0.014 0.018

0.000 0.000 0.000

0.014 0.014 0.017

0.356 0.376 0.262

0.057 0.058 0.071

0.954 0.962 0.951

0.051 0.052 0.057

Table 2 Performance statistics of the slope estimators under Poisson distributions and outlier contamination. Abbreviations: APC, average percent change; RAE, relative absolute error; SD, standard deviation; MSE, mean squared error; ASE, average standard error; AL, confidence interval average length, RR, rejection rate; LM, linear regression; GLM, Poisson regression; QR, median regression. RAE

SD

MSE

ASE

Power

AL

Coverage

RR

Steep/small

LM GLM QR

5.228 4.574 5.305

0.582 0.503 0.604

0.339 0.253 0.364

0.544 0.486 0.598

0.054 0.058 0.048

2.664 2.380 2.928

0.951 0.983 0.952

0.049 0.051 0.047

Steep/large

LM GLM QR

0.781 0.754 0.312

0.107 0.101 0.038

0.011 0.010 0.001

0.080 0.024 0.078

0.456 0.852 0.379

0.392 0.117 0.383

0.968 0.652 0.972

0.036 0.395 0.032

Medium/medium

LM GLM QR

1.786 1.723 1.547

0.107 0.102 0.091

0.011 0.010 0.008

0.099 0.066 0.107

0.076 0.232 0.069

0.429 0.287 0.461

0.956 0.848 0.960

0.045 0.174 0.037

Medium/large

LM GLM QR

1.282 1.224 0.456

0.085 0.079 0.027

0.007 0.006 0.001

0.070 0.019 0.039

0.235 0.713 0.336

0.302 0.082 0.169

0.960 0.534 0.961

0.031 0.496 0.039

Shallow/large

LM GLM QR

1.979 1.865 0.675

0.061 0.057 0.020

0.004 0.003 0.000

0.056 0.014 0.020

0.097 0.622 0.228

0.229 0.058 0.082

0.968 0.447 0.946

0.034 0.567 0.049

APC/mean

Inline Supplementary Table S1 can be found online at http:// dx.doi.org/10.1016/j.canep.2013.08.002. Inline Supplementary Table S2 can be found online at http:// dx.doi.org/10.1016/j.canep.2013.08.002. Table 1 contrasts LM, GLM and QR for different APC and Poisson mean values. For the uncontaminated Poisson distribution, all estimators showed a rejection rate close to the nominal 5%, with

the exception of QR which had a slightly conservative rate at T = 15. GLM outperformed all the other models in terms of RAE, SD and thus MSE for small values of the Poisson mean. For increasing average numbers of events, LM and GLM showed similar results. The median regression, as expected, was slightly less efficient than LM and GLM. However, the MSE was virtually zero for all models at large values of the mean, regardless of the size of the APC. QR

Table 3 Outlier contamination. Deciles of the estimates of the relative change for different time periods at the largest simulated value of the mean parameter. The expected overall relative change is 2. Abbreviations: LM, linear regression; GLM, Poisson regression; QR, median regression. 8 years

Min 10th 20th 30th 40th 50th 60th 70th 80th 90th Max

15 years

30 years

LM

GLM

QR

LM

GLM

QR

LM

GLM

QR

0.02 0.80 1.37 1.63 1.82 1.98 2.17 2.41 2.88 5.10 318.61

0.03 0.84 1.37 1.64 1.82 1.98 2.17 2.42 2.89 5.08 235.98

0.10 1.43 1.62 1.75 1.87 1.99 2.11 2.26 2.46 2.77 49.82

0.00 0.45 0.92 1.36 1.72 2.03 2.42 3.01 4.49 8.99 549.19

0.01 0.47 0.94 1.37 1.72 2.03 2.42 3.00 4.38 8.71 283.23

0.12 1.20 1.46 1.65 1.83 2.00 2.21 2.46 2.78 3.32 11.55

0.00 0.22 0.54 0.94 1.40 2.01 2.96 4.42 7.93 19.87 2333.92

0.00 0.25 0.58 0.98 1.43 1.99 2.87 4.42 7.63 17.23 2282.89

0.19 0.95 1.24 1.49 1.75 2.01 2.34 2.73 3.30 4.28 14.91

M. Geraci et al. / Cancer Epidemiology 37 (2013) 843–849

847

Fig. 2. English cancer registry data, 1990–2006. Time plots of the log-rates of selected tumors for different age groups (first five rows) with regression lines superimposed (linear model, dashed; log-linear model, dotted; median regression model, solid) and density of the log-rates (bottom row) for each age group (darker to lighter shades of grey for, respectively, younger to older groups). Abbreviations: AC, astrocytomas; ALL, acute lymphoid leukaemia; AML, acute myeloid leukaemia; ODG, oligodendroglioma.

848

M. Geraci et al. / Cancer Epidemiology 37 (2013) 843–849

showed a minor overestimation of the standard error at medium values of the Poisson mean being the ASE higher than expected. As a consequence, the AL of the 95% confidence interval was larger than in the other models. All models had coverage close to the nominal 95%, although GLM showed a consistently higher frequency. As expected the situation overturns in favor of the median regression when outliers are introduced (Table 2). The advantage of employing a robust estimator of the slope of the regression line becomes more apparent at higher values of the mean parameter. Both the RAE and the MSE favored QR. The seeming consistently higher power of GLM is, in fact, a consequence of an underestimation of the standard error. It follows that the associated test statistics becomes extremely liberal and that confidence intervals tend to be narrower than expected, as confirmed by the small AL and poor coverage. On the other hand, LM has an error rate comparable to that of QR across all simulation settings. However, the LS estimator of the log-rates showed a higher MSE and a lower power as compared to the LAD estimator at larger mean values. Of course, these results depend on how outlying values were generated. Table 3 shows the impact of outlier contamination in more depth. The cumulative distribution of the APC and MPC estimates was calculated for the entire time periods (8, 15 and 30 years) and compared to the expected doubling of the rates. All models hit the mark at the middle of the distribution. However, the estimators of the average change (LM, GLM) showed a heavy-tailed distribution, thus implying a substantial frequency of APC values away from the middle of the distribution. On the contrary, the distribution of the MPC was denser around its center. For instance, in this simulated scenario, there would be a 30% probability of overestimating the percent change over a 15-year time period by a factor of 150% (or higher) when using the mean estimator. This probability would halve when using the median estimator. Similarly, the probability of underestimating the percent change of a given amount is higher for the mean than for the median regression. These differences became more or less pronounced depending on the span of the time period and, ultimately, on the size of the percent change. Inline Supplementary Table S3 shows performance statistics of the slope estimators under overdispersed distributions. In this case, all models had MSE approximately equal to zero at larger values of the mean. GLM seemed to have an advantage over NB in terms of power. However, the comparison here is complicated by their inflated type I error rates, far from the nominal 5%. Whereas the RR for NB decreased for increasing number of years, GLM showed no improvement. On the other hand, the tests associated with LM and QR had lower power but more reasonable RRs. Inline Supplementary Table S3 can be found online at http:// dx.doi.org/10.1016/j.canep.2013.08.002. 4. English cancer registry data Secular trends in cancer incidence can provide insightful and important clues to the understanding of the etiology of cancer, useful quantitative measures of the effectiveness of campaigns on prevention, as well as a computational base for predicting the future load on national health systems. Statistical analyses often focus on specific age groups for which targeted healthcare services must be provided. For example, this is the case for children or young adults affected by cancer. Although tumors in young people represent a major source of morbidity and mortality [21,22], specific cancer types often pose a challenge in statistical terms due to their low incidence among the population, especially when stratification by age and sex is introduced in the analysis. It is not uncommon in these situations to observe zero counts and, more worryingly, values that depart substantially from the overall pattern. Whether the result of natural variability or of registration practices, such outlying observations might impact on trend estimation.

To illustrate the methods described in the previous sections, we analyzed data provided by the Northern and Yorkshire Cancer Registry and Information Service on cancer incidence in England, 1990–2006. We selected cases aged less than 40 years and assigned them to five age groups. Tumours were classified according to a morphology-based diagnostic scheme [23]. For illustrative purposes, we report the analysis for selected haematopoietic and brain tumours: acute lymphoid leukaemia (ALL), acute myeloid leukaemia (AML), astrocytomas (AC) and ODG. Population estimates were obtained from the Office for National Statistics and yearly rates were calculated per hundred thousand person-years. We then estimated the average and median percent change. In view of the simulation results, the standard error of the MPC was calculated using the ‘rank’ method for all groups except ODG cases aged 20–24 years, for which the ‘nid’ method was employed. Fig. 2 shows the time plots and the estimated density of the logrates by age and diagnostic group, while APC and MPC estimates and their 95% confidence intervals are reported in Inline Supplementary Table S4. Except in very few cases, the density of the transformed rates can hardly be approximated by a normal curve. For ALL cases aged 25–29 years and AML cases aged 20–24 years, the regression lines estimated by LM and GLM crossed the line estimated by QR. As a consequence, the APC and MPC had opposite signs (Inline Supplementary Table S4). The confidence intervals, however, were too wide to conclude whether there was an increase or a decrease of ALL and AML rates in those age groups. In contrast, the MPC for AML cases aged 25–29 years was equal to 0.9% (0.5–2%), whereas the APC had a similar value but wider confidence intervals (Inline Supplementary Table S4). Inline Supplementary Table S4 can be found online at http:// dx.doi.org/10.1016/j.canep.2013.08.002. The number of AC cases in young people aged less than 15 years jumped from approximately 85 cases in the first two years to 120 in the following year, resulting in a rather unusual pattern (Fig. 2). AC rates in this age group increased significantly according to the APC estimates. The latter were equal to 1.8% (0.6–3.1%) and to 1.7% (0.7–2.6%) for LM and GLM, respectively. However, the MPC was not significant and equal to 0.8% (0.0–3%). In this case, the two points at the beginning of the time period have a substantial leverage on the percent change estimated by LM or GLM but they have little effect on the MPC. This should prompt further investigation into possible underlying explanations, such as, for example, changes in coding practice. As reported previously in Fig. 1, ODG rates in people aged 25–29 years showed a high leverage point at the beginning of the time period, where the outliers usually have more weight on the slope of regression lines. After excluding the year 1990 from the calculation, the estimate of the APC was lower (around 2%) and statistically non-significant. 5. Conclusion We presented a novel approach to the estimation of the percent change in incidence and mortality rates over time. The median has robustness properties as well as mathematical properties which have useful practical implications. In an extensive simulation study, our method showed a superior performance as compared to those methods based on LS estimation when the data are contaminated by outliers. In other scenarios, the MPC was less competitive than APC although models that were expected to perform best were affected by the size of the annual rate of change and by the average number of events. We showed that it is difficult to support normality assumptions when few time points are available and the number of events is small as in the case of the cancer data analysis. Another advantage of our method is that median regression does not depend on distributional assumptions.

M. Geraci et al. / Cancer Epidemiology 37 (2013) 843–849

The simulation study provided a useful base of evaluation for APC estimation, in scenarios that can be commonly found in real data applications. For increasing number of events and number of time periods, different asymptotic approximations come into play. Whether the percent change is affected by high leverage points needs to be assessed case by case. Likewise, the choice of calculation method of the MPC standard errors needs to be made accordingly. We believe that, as a good practice, both APC and MPC should be reported when sensitivity issues arise. As John Tukey, one of the founders of modern statistical methods, once commented on the use of robust methods, ‘‘which robust/resistant methods you use is not important – what is important is that you use some. It is perfectly proper to use both classical and robust/resistant methods routinely, and only worry when they differ enough to matter. But when they differ, you should think hard’’ [24, p. 1645]. For example, based on APC results alone, one would conclude that between 1990 and 2006 there was a statistically significant increase of about 1.8% in AC rates among English children and teenagers. In contrast, not only the MPC estimate was much smaller (0.8%) but also non-significant, thus leading to a rather different conclusion. While visual inspection of the data might have a diagnostic value, it does not provide any quantitative support. The MPC is the robust equivalent of the APC which makes its interpretation straightforward. Its calculation, as shown in the supplementary material available with this paper, can be easily carried out with commonly available statistical software, requiring very little effort in terms of programming. Finally, the present study was aimed at addressing the violation of some distributional assumptions under the hypothesis of constant change. Model misspecification of the linear predictor represents a different, although important, issue. For example, in case of changes in trends one could use a spline-type approach [1]. Piecewise-linear robust regression may well offer a possible development for the methods proposed in this paper. Conflict of interest statement None declared. Acknowledgements The Centre for Paediatric Epidemiology and Biostatistics benefits from funding support from the Medical Research Council in its capacity as the MRC Centre of Epidemiology for Child Health (G0400546). The UCL Institute of Child Health receives a

849

proportion of funding from the Department of Healths NIHR Biomedical Research Centres funding scheme. This work was also partially supported by Cancer Research UK. Jillian Birch is Cancer Research UK Professorial Fellow, University of Manchester. We thank three anonymous reviewers for their helpful comments on an earlier version of the article. References [1] Kim HJ, Fay MP, Feuer EJ, Midthune DN. Permutation tests for joinpoint regression with applications to cancer rates. Stat Med 2000;19(3):335–51. [2] Fay MP, Tiwari RC, Feuer EJ, Zou Z. Estimating average annual percent change for disease rates without assuming constant change. Biometrics 2006;62(3):847–54. [3] O’Hara RB, Kotze DJ. Do not log-transform count data. Methods Ecol Evol 2010;1(2):118–22. [4] Huber PJ. Robust statistics. 1st ed. New York, NY: John Wiley & Sons, 1981. [5] Koenker R, Bassett G. Regression quantiles. Econometrica 1978;46(1):33–50. [6] Koenker R. Quantile regression. New York, NY: Cambridge University Press, 2005. [7] R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Core Team, 2013. [8] Koenker R. quantreg: quantile regression; 2012, R package version 4.91. [9] Geraci M. lqmm: Linear quantile mixed models; 2013, R package version 1.03. [10] Anscombe FJ. The transformation of Poisson, binomial and negative-binomial data. Biometrika 1948;35(3/4):246–54. [11] Berry DA. Logarithmic transformations in ANOVA. Biometrics 1987;43(2):439–56. [12] Stewart-Oaten A. Sequential estimation of log(abundance). Biometrics 1996;52(1):38–49. [13] Machado JAF, Santos Silva JMC. Quantiles for counts. J Am Stat Assoc 2005;100(472):1226–37. [14] Powell JL. Least absolute deviations estimation for the censored regression model. J Econometr 1984;25(3):303–25. [15] Powell JL. Censored regression quantiles. J Econometr 1986;32(1):143–55. [16] Powell J. Estimation of monotonic regression models under quantile restrictions. In: Engle R, McFadden D, (, eds. Nonparametric and semiparametric methods in econometrics. New York, NY: North-Holland, 1991 : 357–84. [17] Koenker R. Confidence intervals for regression quantiles. In: Mandl P, Huskova M, eds. Asymptotic statistics. New York, NY: Springer-Verlag, 1994: 349–59. [18] Buchinsky M. Estimating the asymptotic covariance matrix for quantile regression models. A Monte Carlo study. J Econometr 1995;68(2):303–38. [19] Redden DT, Fernandez JR, Allison DB. A simple significance test for quantile regression. Stat Med 2004;23(16):2587–97. [20] Kocherginsky M, He X, Mu Y. Practical confidence intervals for regression quantiles. J Comput Graph Stat 2005;14(1):41–55. [21] Geraci M, Birch JM, Alston RD, Moran A, Eden TOB. Cancer mortality in 13 to 29-year-olds in England and Wales, 1981–2005. Br J Cancer 2007;97(11):1588–94. [22] Alston RD, Geraci M, Eden TOB, Moran A, Rowan S, Birch JM. Changes in cancer incidence in teenagers and young adults (ages 13 to 24 years) in England 1979–2003. Cancer 2008;113(10):2807–15. [23] Birch JM, Alston RD, Kelsey AM, Quinn MJ, Babb P, McNally RJQ. Classification and incidence of cancers in adolescents and young adults in England 1979– 1997. Br J Cancer 2002;87(11):1267–74. [24] Huber PJ, John W. Tukey’s contributions to robust statistics. Ann Stat 2002;30(6):1640–8.