Smoothed Lexis Diagrams With Applications to Lung and Breast Cancer Trends in Taiwan

Cancer surveillance research often begins with a rate matrix, also called a Lexis diagram, of cancer incidence derived from cancer registry and census data. Lexis diagrams with 3- or 5-year intervals for age group and for calendar year of diagnosis are often considered. This simple smoothing approach suffers from a significant limitation; important details useful in studying time trends may be lost in the averaging process involved in generating a summary rate. This article constructs a smoothed Lexis diagram and indicates its use in cancer surveillance research. Specifically, we use a Poisson model to describe the relationship between the number of new cases, the number of people at risk, and a smoothly varying incidence rate for the study of the incidence rate function. Based on the Poisson model, we use the standard Lexis diagram to introduce priors through the coefficients of Bernstein polynomials and propose a Bayesian approach to construct a smoothed Lexis diagram for the study of the effects of age, period, and cohort on incidence rates in terms of straightforward graphical displays. These include the age-specific rates by year of birth, age-specific rates by year of diagnosis, year-specific rates by age of diagnosis, and cohort-specific rates by age of diagnosis. We illustrate our approach by studying the trends in lung and breast cancer incidence in Taiwan. We find that for nearly every age group the incidence rates for lung adenocarcinoma and female invasive breast cancer increased rapidly in the past two decades and those for male lung squamous cell carcinoma started to decrease, which is consistent with the decline in the male smoking rate that began in 1985. Since the analyses indicate strong age, period, and cohort effects, it seems that both lung cancer and breast cancer will become more important public health problems in Taiwan. Supplementary materials for this article are available online.


INTRODUCTION
Cancer surveillance research often begins with a rate matrix, also called a Lexis diagram, of cancer incidence derived from cancer registry and census data. Given an age group and a group of calendar year of diagnosis, the Lexis diagram reports the incidence rate of a disease in terms of the number of new cases per 100,000 person-years. Lexis diagrams with 5-year intervals for age group and 5-year intervals for calendar year of diagnosis are often considered so as to make the data less volatile. Various graphical and quantitative methods have been developed to analyze these Lexis diagrams to study the effects of age, period, and cohort on cancer incidence rates (Anderson 2009;Rosenberg and Anderson 2011). These include classical descriptive plots (Keiding 1990;Devesa, Donaldson, and Fears 1995;Carstensen 2007), age-standardized rates (Last, Abramson, and International Epidemiological Association 1995), estimated annual percentage change (Fay et al. 2006), the joinpoint regression model (Kim et al. 2000), and age-period-cohort models (Clayton and Schifflers 1987a, b;Holford 1991). An excellent example of the use of these methods concluded that the decline Li-Chu Chien and Yuh-Jenn Wu contributed equally; both are co-first authors. Li-Chu Chien, Division of Biostatistics and Bioinformatics, Institute of Population Health Sciences,National Health Research Institutes,. Yuh-Jenn Wu is Assistant Professor, Department of Applied Mathematics, Chung Yuan Christian University, Taiwan (E-mail: yuh-jenn@cycu.edu.tw). Chao A. Hsiung is Distinguished Investigator and Director, Institute of Population Health Sciences, National Health Research Institutes, Taiwan (E-mail: hsiung@nhri.org.tw). Lu-Hai Wang is Distinguished Investigator and Director, Institute of Molecular and Genomic Medicine, National Health Research Institutes, Taiwan (E-mail: lu-hai.wang@nhri.org.tw). I-Shou Chang is Investigator, National Institute of Cancer Research, National Health Research Institutes, Taiwan (E-mail: ischang@nhri.org.tw). Dr. Lu-Hai Wang's work is supported by the Center of Excellence in Cancer Research, Department of Health, Taiwan (DOH102-TD-C-111-004). Dr. I-Shou Chang's work is partially supported by grants from the National Science Council (NSC103-2811-M-400-005) and the Ministry of Health and Welfare (MOHW103-TDU-212-114001). The authors thank three referees, an associate editor, and the editor for their critical reading and constructive suggestions that improve the article. of lung cancer incidence around 2000 among U.S. women is likely to persist through at least 2025, based on calendar year and birth-cohort age-specific incidence patterns (Jemal, Ward, and Thun 2005).
While these methods are useful in generating hypotheses and have contributed significantly to cancer studies, they often start with a 5-year, instead of 1-year, tabulated Lexis diagram. Using 5-year tabulation represents one intuitive smoothing method to reduce the volatility in the data and to facilitate subsequent studies. This simple smoothing approach suffers from a significant limitation; important details may be lost in the averaging process involved in generating a summary rate, and these details may be useful in understanding time trends in disease (Holford 1991).
The purpose of this article is to present a smoothed Lexis diagram and use it to study the effects of age, period, and cohort on incidence rates in terms of straightforward graphical displays. These include the age-specific rates by year of birth, age-specific rates by year of diagnosis, year-specific rates by age of diagnosis, and cohort-specific rates by age of diagnosis (Carstensen 2007). We will illustrate our approach by studying the trends in lung and breast cancer incidence in Taiwan.
Our approach capitalizes on the assumption that cancer incidence rate varies smoothly with both age and calendar time. Cancer arises from combinations of changes that occur in the same cell over a long period of time. Only when a critical number of changes occur in the same cell does it finally become cancerous. Known causes of cancer include genetic factors; lifestyle factors such as tobacco use, diet, and physical activity; and environmental exposure to different types of chemicals and radiation. The above assumption amounts to saying that causes of cancer, known or unknown, do not change abruptly with age or calendar year. This also has something to do with cancer screening projects, which we comment on in the discussion.
To simplify the discussion, we assume that the unit of time is a year, which is a positive integer, and that we consider a specific cancer. Let M(x, y) denote the number of people not having been diagnosed with this cancer before age x and year y, and N(x, y) denote the number of people who are newly diagnosed with this cancer at age x and year y. With proper transformation, we assume both x and y are represented by fraction numbers in [0,1]. Let F (x, y) denote the probability that an individual at age x will be a newly diagnosed case in year y. F (x, y) will be referred to as incidence rate. The main assumption of this article is that N(x, y) is a Poisson random variable having parameter M(x, y)F (x, y); thus, We note that when the incidence rate F (x, y) is small and the risk set size M(x, y) is large, the Poisson model approximates the binomial model with parameters M(x, y) and F (x, y); this is the well-known law of rare events (Taylor and Karlin 1984). We also note that, while the binomial model is more intuitive, the Poisson model is computationally more favorable when M(x, y) is large and F (x, y) is small.
Assuming that the incidence rate F varies smoothly with year of diagnosis and age of diagnosis, our main task is to estimate F based on M(x k , y k ) and N(x k , y k ) for k = 1, . . . , K. In applications, (x k , y k ) are often lattice points in [0, 1] × [0, 1], M(x k , y k ) and N(x k , y k ) can be obtained from data in census bureau and cancer registry, and K is the total number of lattice points considered.
Since F is a bounded nonnegative function, this is a shaperestricted inference problem. We will take a Bayesian approach that puts a prior on the space of bounded nonnegative smooth functions on [0, 1] × [0, 1]. The age x-specific incidence rate by calendar year y of diagnosis will be reported based on the posterior distribution of F (x, y). The trends of the incidence rate can be studied by the estimates of F (x, y). We note that since Bernstein polynomials are differentiable, we can even explore the trends of the incidence rate by their partial derivatives, which often provide additional insights.
The prior will be specified in terms of two-dimensional Bernstein polynomials. We note that one-dimensional Bernstein polynomials have been successfully used in Bayesian shaperestricted inference problems regarding a smooth function on [0, 1], using the facts that priors on Bernstein polynomials can be introduced through their coefficients, the geometry of a Bernstein polynomial can be read off readily from its coefficients, and the support of the space of Bernstein polynomials can be large (Chang et al. , 2007Chien et al. 2009). We now indicate that these properties can be extended to two-dimensional Bernstein polynomials and can be used to introduce priors for Bayesian shape-restricted inference for smooth functions on the unit square.
This article is organized as follows. Section 2 presents the statistical model using two-dimensional Bernstein polynomials, proposes the prior distribution, and describes the main observations used in the algorithm for sampling the posterior distribution. The algorithm itself is given in Appendix A. Section 3 presents a simulation study, which mimics one of the real datasets and paves the way for the real data analysis in Section 4. To illustrate our method, we study the trends of lung and breast cancers in Taiwan in Section 4. More detailed information about these cancers in Taiwan is given in Appendixes B-F. Section 5 provides a brief discussion. All Appendixes are given in the supplementary materials.
With Proposition 1, Proposition 2, and (1), we are motivated to introduce a prior distribution on I by defining a probability distribution on B. Let π be a probability measure on B. π is called a Bernstein prior. A natural way to introduce a prior π is to consider π (n, b n ) = π n (b n )p(n), with ∞ n=0 p(n) = 1 and π n a density function on B n ; in fact, π n (•) is the conditional density of π on B n and also denoted by π (•|{n} × B n ). The following proposition shows that the support of the Bernstein priors can be quite large.
Proposition 3. LetB n = {b n ∈ B n : F b n ∈ I n }. Assume p(n) > 0 and the conditional density π (b n |{n} ×B n ) has sup-portB n for infinitely many n. Let F be a continuous realvalued function defined on [0, 1] The proof for Proposition 3 is similar to that for the Proposition 3 in Chang et al. (2007) and hence omitted.
The second observation is that if F b n (x, y) is close to the true incidence F (x, y), then the raw incidence rates (Lexis diagram) N (x k , y k )/M(x k , y k ) for k = 1, . . . , K provide some idea of the range of the values of F b n . This, in turn, can be used to propose priors on the coefficients, b n = (b 0,0,n , . . . , b n,n,n ), in view of Proposition 1.

Sampling the Posterior Distributions
Since the parameter space consists of subspaces of different dimension, we will use a Metropolis-Hastings reversible jump (MHRJ) algorithm (Green 1995) to generate data from the posterior distributions. While we will postpone the algorithm to Appendix A, which is technical and long, we will present here one observation that motivates the way we update the coefficients of the Bernstein polynomials. Since the state space of our Markov chain Monte Carlo (MCMC) algorithm consists of Bernstein polynomials, the closeness concept suggested by the Bernstein polynomial approximation theorem is useful in de- , even for m = n. This suggests that in updating the Markov chain, it is a good idea to make the new Bernstein polynomial "close" to F b n (•, •) even if m = n. There are several obvious methods to use this closeness concept in the algorithm. The one in our algorithm seems to be more computationally efficient, according to our preliminary simulation studies (not reported).
The prior in this simulation study and the following real data analyses assume p(0) = p(1) = 0, which means only Bernstein polynomials of order 2 or higher are included. This is because our preliminary analyses showed that the results obtained under the assumption p(0) = p(1) = 0 are all very similar to those under p(0) = 0, but the Markov chain converges much faster under p(0) = p(1) = 0.
We use the MHRJ algorithm with c = 0.35 to generate the posterior distribution. This value, c, controls the probabilities of changing the order of the polynomial in updating the chain; and c = 0.35 allows a fairly large such probability.
In this article, we take the following strategy to sample from the posterior distribution. Given a set of priors, we run five MCMC with initial values chosen randomly from the priors and monitor convergence by the Gelman-Rubin statisticR, suggested in Gelman and Rubin (1992) and Gelman et al. (2004). The Gelman-Rubin statisticR is calculated for the volume under the surface. We first run each chain with 200,000 iterations and use the second halves of the five chains to computeR. If R is smaller than 1.1, we stop. Otherwise, we double all chain lengths and again use the second halves to computeR. We repeat the procedure untilR is less than 1.1; the updates from the second halves of these five sequences are considered sample from the posterior distribution, which form the basis for inference.
The estimate of the incidence rate function is defined to be the posterior mean and is denoted byF . One naïve way to obtainF is to compute, for each of the 1020 grids, the mean of the values at the grid of the Bernstein polynomials in the updates. In case the orders of the polynomials in the posterior distribution are constant,F is simply the mean of these Bernstein polynomials. We note that in all the Bayesian analyses in this article the polynomial order in the posterior distribution is no larger than 3; in most of the analyses, the polynomial order is exactly 2. Table 1 describes certain properties of the estimates in the simulation studies. It indicates that the numerical performance of the Bayesian method in this study is quite good. Each of the rows presents the performance under one set of priors. Columns 2 and 3 give, respectively, the L 1 -norm and sup-norm of |F − F T | × 10 5 , where | • | represents the absolute value. Columns 4 and 5 give, respectively, the posterior mean and maximum ofF × 10 5 . Columns 6, 7, and 8 report, respectively, the Gelman-Rubin statisticR, the number of iterations that achieve it, and the CPU time for this computation. Column 9 reports the posterior predictive p-value, where the discrepancy function is the maximum of the difference between the observed incidence rates and the incidence rates specified by the model. There are other reasonable discrepancy functions to use, for example, we may use the one that measures the average of the difference between the observed rates and the rates specified by the model. We note that a large posterior predictive p-value suggests no inconsistency between the model and the dataset; see Gelman et al. (2004). Table 1 shows that our method provides a good estimate of the true incidence rates, is insensitive to the choice of prior distributions, and is computationally efficient. In fact, using the criterion T = max{|F priori (x k , y k ) −F priorj (x k , y k )||k = 1, . . . , 1020} × 10 5 , the largest difference happens between set of priors 1 and 4 and has the value of 4.711. Figure 1 shows the simulation study results with F T given in (2). The upper panel shows the year-specific incidence rates by age of diagnosis (left); age-specific incidence rates by year of diagnosis (right). The lower panel shows the absolute value of the difference of year-specific incidence rates by age of diagnosis between the estimateF and the true incidence rate F T multiplied by 10 5 (left); that of age-specific incidence rates by year of diagnosis between the estimateF and the true incidence rate F T multiplied by 10 5 (right).
We now use this simulation study to compare our Bayesian method to the standard Lexis diagram in estimating the true incidence rate F T . The 1-year tabulated standard Lexis diagram consists of the r k = z k M(x k , y k ), for k = 1, . . . , 1020, described above. Based on the 1-year tabulated standard Lexis diagram, we construct 3-, 5-, and 7-year tabulated standard Lexis diagrams in the same manner as that described in Carstensen (2007); see Table 1 and Figure 2 in Carstensen (2007). We may shorten the name of these Lexis diagrams, for example, the 5-year tabulated standard Lexis diagram may be referred to as the 5-year Lexis diagram. These standard Lexis diagrams provide estimatesF of the true incidence rate F T ; it is reasonable to regard them as estimates at the middle points of the Lexis diagram. For the 1-year standard Lexis diagram, the middle points are defined to be the 1020 lattice points mentioned above. For other standard Lexis diagrams, the coordinates of the middle points are defined to be the averages of the coordinates of the lattice points involved in the averaging process that constructs the standard Lexis diagram. (This definition of middle points is slightly different from that in Carstensen (2007), but is more appropriate here in view of our data generation procedure.) Thus, when comparing our Bayesian method to a standard Lexis diagram, we consider the L 1 -norm and sup-norm of |F − F T | × 10 5 over the middle points defined by each respective Lexis diagram; these comparison results are reported in Table 2. For this comparison, we use the set of priors number 4 in Table 1, because it has the largest posterior predictive p-value. Table 2 indicates clearly that our Bayesian method outperforms standard Lexis diagrams in providing estimates of the incidence rates. Our Bayesian method has the additional advantage of producing smooth estimates at every point in the unit square. We note that Figure 1(b) suggests that the largest error in the estimation of F occurs at the boundary.
The results shown in Tables 1 and 2 are quite general. In fact, we carried out several such studies with different F T and M(x, y) and obtained similar results. In the interest of conciseness, we present in Tables 3 and 4 one more study with the same M(x, y) but a different F T defined by In fact, the only difference between these two studies is that the true incidence rate function F T in (3) is not monotone in y. Tables 3 and 4 are, respectively, the counterparts of Tables 1 and 2. Together, Tables 3 and 4 show again that our Bayesian method outperforms standard Lexis diagrams.

Preliminaries
Taiwan is a densely populated island with one of the fastest growing aging populations in the world and has experienced tremendous lifestyle changes in the past several decades. With cancer rates on the rise, it is of interest to investigate both the time trends and secular trends to monitor changes in the incidence rates of various cancers. For example, using the Taiwan Cancer Registry from 1980 to 2006, in which age-standardized incidence rates of most important cancers were reported, the effects of the 1995 implementation of Taiwan's National Health Insurance on cancer incidence rates were explored (Chiang et al. 2010).
There also appeared more detailed epidemiological studies of specific cancers. For example, using the Taiwan Cancer Registry and data from SEER program, the incidence of FIBC (female invasive breast cancer) in the younger generation of Taiwanese was found to be rapidly increasing and approaching that of Caucasian Americans, suggesting a consequence of Westernization of lifestyles (Shen et al. 2005). As another example, lung ADC (adenocarcinoma) was found to be the predominant lung cancer histology subtype among Taiwan female patients; furthermore, its incidence rate was increasing, based on data from National Taiwan University Hospital (Chen et al. 2005).
As an illustration of the method of this article, we will study age, period, and cohort effects on breast cancer and lung cancer in Taiwan. Specifically, we consider people with age between 35 and 85 years in the period from 1988 to 2007, based on Figure 1. This figure provides graphs of the estimateF and the absolute difference between the estimateF and the true incidence rate F T . the Taiwan Cancer Registry and Taiwan census data. Since the incidence rate of female lung SCC (squamous cell carcinoma) is very low, we consider only FIBC, female lung ADC, male lung ADC, and male lung SCC.
The Taiwan Cancer Registry was launched in 1979; it collects information on newly diagnosed cancer patients from all hospitals in Taiwan with a capacity of more than 50 beds. Knowing that the quality of the data has improved steadily, we focused the current study on a recent 20-year period. For a more detailed description of the Taiwan Cancer Registry, see Chiang et al. (2010).

Construction of the Smoothed Lexis Diagram
We now describe the way we construct the smoothed Lexis diagram for the incidence of female lung ADC. The smoothed Lexis diagrams of other cancers are constructed similarly. Using the notation of Section 3, for the lattice point (x k , y k ) in [0, 1] × [0, 1], the number of newly diagnosed female lung ADC patients at the age 35 + 50x k in the year 1988 + 19y k is denoted by N (x k , y k ), which were obtained from the Taiwan Cancer Registry on June 11, 2010. M(x k , y k ) is similarly defined and was derived jointly from Taiwan census data and the Taiwan Cancer Registry. The ICD-9-CM codes for the cancer site (topography) and the ICD-O3 codes for the histology (morphology) used in this study are given in Appendix B.
Let r k = N (x k , y k ) M(x k , y k ) for k = 1, . . . , 1020; we analyze the real datasets in the same manner as we conduct simulation studies in Section 3. In particular, we consider the same four sets of priors, carry out the same MCMC algorithm for posterior distribution sampling, use the Gelman-Rubin statistic to monitor the convergence, and apply the posterior predictive p-value to assess the model fitness.
Tables C1-C4 in Appendix C summarize the MCMC computation for FIBC, female lung ADC, male lung ADC, and male lung SCC, respectively. Entries in these tables bear the same meaning as those in Table 1. To conduct sensitivity analyses, we compare the results obtained from analyses using the four different sets of priors for each cancer dataset. Tables C1-1, C2-1, C3-1, and C4-1, also in Appendix C, report the differences between the posterior means studied in Tables C1, C2, C3, and C4, respectively. These tables indicate that, except in the case of male lung SCC, the models considered are quite consistent with the datasets and the results are not sensitive to the priors. In fact, even for male lung SCC, two of the sets of priors show acceptable posterior predictive p-values and the posterior means for these two priors are quite close to each other. The following results are based on those priors showing largest posterior predictive p-values. Figure 2 presents the smoothed Lexis diagrams, with 2a, 2b, 2c, and 2d for FIBC, female lung ADC, male lung ADC, and male lung SCC, respectively. Contour lines with associated incidence rates are shown in these figures; by definition, incidence rates at points on the same contour line are the same; they are given in terms of number of cases per 100,000 person-years. For example, Figure 2(c) shows that in the year 1988, the incidence rate for an 83-year-old male was about 60 cases per 100,000 person-years. These contour lines help in grasping the pattern change of the incidence rate for each cancer. For example, while the incidence rate of lung ADC increases with both age and calendar year, male lung SCC incidence rate peaks around 2003 among older people. In view of the reports that smoking is a risk factor for lung cancer, more for SCC than for ADC, and that among never smoking lung cancer patients, ADC subtype is more prevalent than SCC (Sun, Schiller, and Gazdar 2007;Couraud et al. 2012), it is of interest to see if the pattern of contour lines in Figure 2 can be explained by a decrease in the smoking rate in Taiwan.

Age, Period, and Cohort Effects on Incidence Rate
We now explore the age, period, and cohort effects on incidence rate in terms of the four classical plots derived from a Lexis diagram (Carstensen 2007). For each of the smoothed Lexis diagrams constructed in last subsection, these four plots are similarly constructed and reported in Figures 3-6. They are, respectively, the plots for age-specific rates by year of birth (rates vs. age of diagnosis, observations within each birth-cohort are connected), age-specific rates by year of diagnosis (rates vs. age of diagnosis, observations within each year of diagnosis are connected), year-specific rates by age of diagnosis (rates vs. year of diagnosis, observations within each age group are connected), and cohort-specific rates by age of diagnosis (rates vs. birth-cohort, observations within each age group are connected).
Each of these figures has four parts, indexed by a, b, c, and d, for FIBC, female lung ADC, male lung ADC, and male lung SCC, respectively.
Since our main purpose is to illustrate the method, we explore some possible uses of each of these figures, even if the same conclusion may be observed from different figures.
The following observations are clear from Figure 3. Given age of diagnosis, later birth-cohorts experience a higher incidence rate for lung ADC. Given a birth-cohort and regarding incidence rate as a function of age, this function is increasing and convex for both female and male lung ADC; it is increasing for male lung SCC. It is also increasing for the younger cohort for FIBC; but it is concave with a peak for birth-cohorts around 1928 and 1923; we expect to see peaks in later cohorts if we observe the rates in longer periods. It might be of interest to see if these observations hold in other geographical regions and to provide an interpretation if they do. Figure 3(b)-(d) indicates that our results are in line with the findings that, among different lung cancer histology types, adenocarcinoma is the most prevalent and its incidence is increasing in East Asia and West Europe (Janssen-Heijnen and Coebergh 2003; Chen et al. 2005;Devesa et al. 2005;Liam et al. 2006;Toyoda et al. 2008). The facts that incidence rate is increasing and convex suggest a rapid increase tendency. One way to stress this is to compute their derivatives. Figure D1 in Appendix D shows that these derivatives for lung ADC are not only positive but also increasing by themselves, and that those for male lung SCC are positive but not always increasing. Figures D1-D4 in Appendix D are, respectively, the derivatives for those in Figures  3-6. Figure 4 indicates that given year of diagnosis and regarding incidence rate as a function of age, this function is increasing for lung ADC and male lung SCC, whereas it is unimodal for FIBC with a peak around 53 years of age. For FIBC and lung ADC, Figure 4 also indicates that, except for male lung SCC, the age-specific rate increases with year of diagnosis for nearly every age group.   . Year-specific rates by age of diagnosis (rates vs. year of diagnosis, observations within each age group are connected). Figure 6. Cohort-specific rates by age of diagnosis (rates vs. birth-cohort, observations within each age group are connected). Figure 5 indicates that for every age specification, incidence rate as a function of year of diagnosis increases for lung ADC; for some age groups, this does not hold for FIBC and male lung SCC. For example, Figure 5(d) shows that the male lung SCC incidence rate first increases and then decreases for older people. This is consistent with the declining smoking rate in Taiwan. We note that for the period from 1971 to 2010, the adult female smoking rate has been fluctuating around 4% and that for males decreased from 60% in 1985 to 35% in 2010, as described in Appendix E. Figure 5 also shows that, for a given year of diagnosis, older people experience a higher incidence rate of lung cancer. For FIBC, peak incidence occurs around the age of 53, which can also be seen from Figure 4(a) and Figure D2a. In fact, Figure D2a suggests that the age that peaks increases slightly with calendar year; this phenomenon deserves further attention, in view of the remark that in Japan and Taiwan the peak age does not seem to be increasing with calendar year as in Hong Kong and the United States (Shen et al. 2005). We note that the incidence rate of FIBC in California, USA, peaks around 70 (Clarke et al. 2012). Figure 6 shows that given an age specification, people in later cohorts experience higher rates, except for male lung SCC and FIBC. For example, Figure 6(d) shows that for males 75 years of age, the 1928 cohort experiences higher lung SCC rates than nearby cohorts. This is also consistent with the declining smoking rate in Taiwan. For FIBC, among females 85 years of age, the 1918 cohort experiences higher FIBC rates than nearby cohorts, but among the younger cohorts, later ones experience higher rates.

General Observations.
We illustrate the use of smoothed Lexis diagrams in studying lung and breast cancer trends in Taiwan by presenting the four classical plots that depict the age, period, and cohort effects on incidence rates. Except for male lung SCC, the smoothed Lexis diagrams are based on the Bayesian models that fit the datasets well and the posterior distributions are not sensitive to the priors (which enhances our confidence in the four classical plots). We now compare those plots with ones obtained from standard Lexis diagrams to examine the gains and limitations of our approach.
Here are some general remarks from this comparison. Oneyear tabulated standard Lexis diagrams are volatile, and the four classical plots based on them are typically noisy. In general, the longer the tabulation period of a standard Lexis diagram is, the smoother the corresponding classical plots, the shorter the domains of definition of these classical plots, and the fewer the number of the classical plots that can be obtained. However, our smoothed Lexis diagram does not suffer any such drawbacks.
The main observation of the comparison based on the four cancers in Taiwan is that our methods reveal cancer incidence trends more clearly than standard Lexis diagrams; if the trends do not involve the oldest age group, then they are not in contradiction with the results from standard Lexis diagrams, which may or may not reveal any trends; for trends regarding the oldest age group, there may exist contradictions, for which extra efforts may offer some insights.
We now present two examples to portray the above observations. The first regards age-specific rates by year of diagnosis for breast cancer. The smoothed Lexis diagram shows that, around the year 2006, people around the age of 53 experienced the highest incidence rate among all age groups, while the 3-year standard Lexis diagram shows that, around the same period, people between the ages of 50 and 65 experienced the highest incidence rate; see Figure 4(a) and Figure F2-3a in the supplementary materials. We note that the age range showing the highest incidence rate gets smaller as the tabulation period gets longer and that the 13-year standard Lexis diagram reports that the highest incidence rate is around 55-years of age (data not shown). We also note that the 13-year Lexis diagram for our data is a 2 × 4 matrix.

4.4.2
Choosing the Age Range. The second example regards a possible discrepancy shown in the plots for age-specific rates by year of diagnosis among the oldest male lung SCC patients. The smoothed Lexis diagram indicates that the incidence rate increases with age for each year of diagnosis, whereas the standard Lexis diagram indicates that for certain years of diagnosis, it increases with age only before age 75 (see Figure  4(d) and Figure F2-3d). Thus, we can more confidently claim that the male lung SCC incidence rate increases with age until at least age 75. To gain some insights into this discrepancy, we examine several standard Lexis diagrams simultaneously. Based on the 3-year Lexis diagram, it seems that incidence rates increase over larger age ranges for more recent years of diagnosis than for earlier years of diagnosis. This pattern becomes more prominent for 5-and 7-year period tabulated standard Lexis diagrams, as well as for longer periods; eventually, the incidence rate increases everywhere in its domain for 13-year Lexis diagram (data not shown). These plots show that the agreement between our results and those from standard Lexis diagrams is higher for more recent years of diagnosis.
The above discrepancy seems to be in line with the small predictive p-values shown in Table C4 and mentioned in Section 4.2. Small predictive p-values suggest inconsistency between model and dataset and may serve as warnings that prompt for more detailed studies. We think one possible source for this discrepancy regards the fact that incidence rate estimation error is likely to be biggest for the oldest people in 1988, because Taiwan is a fast-growing aging population with the number of male at age 85 being 4103 for 1988 and being 19,630 for 2007. To avoid this type of difficulties, we considered smoothed Lexis diagrams for people with different age ranges, including [35,65], [35,70], [35,75], [65,85], and [70,85]. We find the posterior predictive p-values are not small and the classical plots based on them are more satisfactory (details are omitted).

Remarks on the Taiwan Lung and Breast Cancer Trends
The single most prominent finding in the above smoothed Lexis diagrams is that disease burden, as a function of year of diagnosis, is steadily increasing in lung ADC and FIBC, and decreasing in male lung SCC, for the period from 1988 to 2007. While the latter is consistent with the reduction of smoking rate, the former deserves serious attention from the point of view of public health.
Although the incidence rates of both lung ADC and FIBC are on the rise, they show different patterns. Figure 4 shows that considering incidence rate as a function of age of diagnosis, given year of diagnosis, lung ADC, and FIBC show very different patterns, with the former increasing in a much larger age range. In our study, FIBC peaks between 50 and 55 years of age; but in the United States it peaks at about age 70 (Shen et al. 2005), which prompts questions on the role of environmental contributions to FIBC susceptibility. In view of the discussion of the impact of Westernization of lifestyles on FIBC incidence rates in Taiwan (Shen et al. 2005) and of the main molecular subtypes of FIBC in Taiwan (Lin et al. 2009), we think it of interest to learn why FIBC peaks at such an early age in Taiwan. From Figure 3(a), we know that except for the 1923-1928 birth cohorts, age-specific incidence rate increases with age for each birth cohort in the period from 1988 to 2007; for birth cohort around 1923-1928, it increases until age 75. We think age-specific rates by year of birth are relevant in addressing the impact of Westernization of lifestyles on FIBC incidence rates and we need to obtain them from other periods and countries to better understand this issue.
Given that the smoking rate among males is decreasing, that the smoking rate among females has remained small, and that the incidence rate of lung ADC is rapidly increasing among both males and females, it is of interest to look for other risk factors for the development of lung ADC in Taiwan (Hsiung et al. 2010;Lan et al. 2012;Lo et al. 2013). Although there have been advances in these directions, further studies are needed to address this public health issue.

DISCUSSION
We use a Poisson model to describe the relationship between the number of new cases, the number of people at risk, and a smoothly varying incidence rate for the study of the incidence rate function. Based on the Poisson model, we use the standard Lexis diagram to introduce priors through the coefficients of Bernstein polynomials and propose a Bayesian approach to construct a smoothed Lexis diagram. We outline a fairly standard procedure for implementing the MCMC algorithm. For the four cancers studied in this article, we generally find that the models fit the datasets well, the results do not seem to be sensitive to the priors, and it is a good idea to construct the four classical plots to describe directly the effects of age, period, and birth-cohort on cancer incidence rates. While these classical plots from smoothed Lexis diagrams are revealing, comparing them with their counterparts based on standard Lexis diagrams and considering smoothed Lexis diagrams based on different age range helps to gain more confidence in the results or to explore subtle patterns regarding older people.
The prior distributions consist of two parts: the polynomial order and the polynomial coefficients. As discussed in Section 2, the geometric properties of Bernstein polynomials suggest we introduce the prior on polynomial coefficients through the raw rate matrix. As for the prior on the polynomial order, we used a truncated Poisson distribution out of convenience. For the simulation study and the applications, it turns out that the posterior distribution of the polynomial order is exactly 2 for most of the studies and is never larger than 3, which are much smaller than the polynomial order in the prior. We think it is a good idea to try several priors on the polynomial order to see if the results are sensitive to the choice of the prior.
Our approach assumes a smoothly varying incidence rate on the age range and period range under study. This assumption may not hold when screening rate changes due to factors such as media coverage of a famous person with cancer (Lane, Polednak, and Burg 1989), a publicity campaign to increase screening (Cram et al. 2003), implementation of a specific cancer screening project , or a nationwide health insurance policy (Chiang et al. 2010). Although some of these screening effects may not affect long-term trends and smoothing over them may not cause too much concern (Feuer and Wun 1992), it is advisable to always pay attention to this assumption when applying the method of this article. We note that our Bayesian analysis includes posterior predictive p-values to assess the fit of the model to the data and small predictive p-values may serve as warnings that further investigations are needed, which may be useful in the choice of age range and period range under study.
Taiwan's national health insurance program started in 1995; it is of interest to see if this event might affect the otherwise possibly smoothly varying incidence rate function and consequently the results in this article. One possible way to deal with this issue is to compare smoothed Lexis diagrams of periods before and after 1995 for each cancer. It would also be of interest to compare these patterns across different cancers. While our example on male lung SCC exemplifies the importance of considering Lexis diagrams covering different ranges of age of diagnosis, this discussion suggests the possible use of Lexis diagrams covering different ranges of year of diagnosis.
While this approach allows formal Bayesian inference, our illustrations focus on graphical displays. Other uses of the posterior distribution of incidence rates should be possible. For example, it might be interesting to explore the posterior distribution to propose age, period, and cohort models for closer surveillance studies.

SUPPLEMENTARY MATERIALS
The supplementary materials include Appendices A-F.

NOTES
All the computer programs in this article are written in Matlab. These include the MCMC algorithm for sampling the posterior distributions and the graphical displays of the smoothed Lexis diagrams. The codes are available from the corresponding author I-Shou Chang.