An analysis of Japanese liver cancer mortality data with
Bayesian age–period–cohort models
Wataru Sakamoto ∗Abstract — Age–period–cohort (APC) models have been widely used in the analysis of incidence and mortality data.
Bayesian APC models, in which multivariate Gaussian priors are incorporated on age, period and cohort eﬀects, can evadethe identifiability problem. Inference with integrated nested Laplace approximations (INLA) has recently been a useful tool.
An application of the Bayesian APC models with INLA to Japanese liver cancer mortality data is illustrated, in which asudden change of the cohort eﬀect was revealed.

Keyword: information criteria; integrated nested Laplace approximation; Gaussian Markov random field
1Introduction

Cancer has been the leading cause of death for Japanese people. To analyze cancer incidence andmortality data and to predict mortality rate in the future should be essential to establishing an evidencebased target in cancer control strategy.
Age–period–cohort (APC) models have been widely used in epidemiology, such as in the analysisof incidence and mortality data. However, classical APC models have the identifiability problem (e.g.
Kupper et al., 1985), because of the identity: period = cohort + age, and some approaches to overcome
the identifiability problem have been proposed. One way is to consider linear constraints amongparameters or parsimonious parametrization (e.g. Carstensen, 2007). R packages Epi (Hills et al., 2014)
and apc (Nielsen, 2015) are available for the APC analysis.
Another approach to evade the identifiability problem is to put regularization constraints to thevariation in each component of age, period and cohort. This corresponds to a Bayesian approach inwhich multivariate Gaussian priors are incorporated on age, period and cohort eﬀects. Nakamura(1986) considered a Bayesian APC model, in which successive parameters on age, period and cohort
eﬀects were assumed to change gradually, and smoothing parameters were estimated by an empirical
Bayes method. Bayesian hierarchical approaches to the APC model have also been studied; see Schmidand Held (2007) (and references therein), in which Markov chain Monte Carlo (MCMC) algorithm wasused for inference. Riebler et al. (2012) used integrated nested Laplace approximations (INLA) (Rue etal., 2009) for inference in Bayesian APC models, which enable fast computation and have equivalent
performance with MCMC.
In this paper we apply Bayesian APC models to Japanese liver cancer mortality data, and illustrate
usefulness of the Bayesian APC models and inference with INLA.
2

Bayesian age–period–cohort models

Let yi jk be the number of deaths or incidences and nijk be the populations, observed for the i-th age
group (i = 1, . . . , I), j-th period group ( j = 1, . . . , J) and k-th cohort group (k = 1, . . . , K). The following
identities hold: Agei + Cohortk = Period j . Consider a Poisson modelyijk ∼ Po(nijk λijk ),i = 1, . . . , I; j = 1, . . . , J; k = 1, . . . , K,

(2.1)

where λi jk is the unknown relative risk. The APC model has the following form on the linear predictor
ηi jk :log(λijk ) ≡ ηijk = µ + αi + β j + γk ,(2.2)
∗ Faculty of Environmental Science and Technology, Okayama University, 3-1-1 Tsushima-naka, Kita-ku, Okayama 700-8530,
Japan, E-mail: w-sakamoto@okayama-u.ac.jp

where µ, α j , β j and γk are unknown constants representing mean, age, group and cohort eﬀects,respectively.
In Bayesian analysis of the APC model (2.2), we incorporate priors such as follows on each of theage, period and cohort eﬀects: for example, on age eﬀect {αi },• Random walk prior of first order (RW1): p(αi ) ∝ 1 marginally, andαi − αi−1 ∼ N(0, θ−1α ),i.e.,

αi |αi−1 ∼ N(αi−1 , θ−1α ),i = 2, . . . , I,

which assumes a constant trend over time. The precision parameter θα controls variation.
• Random walk prior of first order (RW2): p(αi ) ∝ 1 marginally, andαi − 2αi−1 + αi−2 ∼ N(0, θ−1α ),i.e.,

αi |αi−1 , αi−2 ∼ N(2αi−1 − αi−2 , θ−1α ),i = 3, . . . , I,

which assumes a linear trend over time. The precision parameter θα controls smoothness.
These priors play a role of solving the identifiability problem (Schmid and Held, 2007). Priors on
period and cohort eﬀects {β j } and {γk } are incorporated similarly, with precision parameters θβ andθγ , respectively. Moreover, hyper-priors such as gamma distributions are introduced to the precisionparameters. For the constant term µ in (2.2), we usually consider a non-informative flat prior.
33.1

Bayesian inference
Latent Gaussian models

The Bayesian APC model described in Section 2 belongs to the class of latent Gaussian models (Rue etal., 2009), which is characterized as follows:• y = {yi jk } is a vector of observations, each of which follows independently the distributionp(yijk |ηi jk , θ) for a given latent variable, ηijk in the APC model, as described next.
• z = {ηijk , µ, αi , β j , γk } is a high-dimensional vector of latent variables, of which the prior forms
an intrinsic Gaussian Markov random field (GMRF) (Rue and Held, 2005), which follows themultivariate normal distribution with mean vector 0 and sparse precision matrix Q(θ). Let p(z|θ)
be its density.
• θ = (θα , θβ , θγ ) is a hyper-parameter vector. Let π(θ) be a hyper-prior on θ, which is usuallynon-normal.
The joint posterior density of (z, θ) becomes, via the Bayes’ theorem,∏π(z, θ|y) ∝p(yijk |ηi jk , θ) · p(z|θ)π(θ).
i,j,k

3.2

Integrated nested Laplace approximation

The latent Gaussian model described above is regarded as a hierarchical Bayes model, and we aremainly interested in predicting each latent variable, say zi , that is, obtaining its marginal posterior
distribution:"p(zi |y) =π(z, θ|y)dz−i dθ.
(3.1)
The marginalization in (3.1) requires high-dimensional integration. The integrated nested Laplaceapproximation (INLA), proposed by Rue et al. (2009), enables fast computation without posterior sampling, using the Laplace approximation of the posterior distribution of the sparse latent variables, which
is normally approximated, and using numerical integration on the distribution of hyper-parameters,
which is non-normal. It has been shown that the INLA has comparable estimation performance to theMarkov chain Monte Carlo method (Rue et al., 2009; Riebler et al., 2012).
An R package INLA (http://www.r-inla.org/) (Martino and Rue, 2010) for Bayesian analysis withthe latent Gaussian models and the INLA is available, and we used it through our study.

4Analysis of Japanese liver cancer mortality data

Japanese cancer mortality data are based on Vital Statistics Japan (Ministry of Health, Labour andWelfare), and can be available from the website of Cancer Registry and Statistics, Cancer Information
Service, National Cancer Center Japan: http://ganjoho.jp/reg_stat/statistics/dl/index.html.
The populations and the numbers of liver cancer death, separated by gender, are recorded every year
from 1958 to 2014 for each of 5-year age categories. The cohort category has 5-year length and is shifted
year-by-year. Analyses of the Japan liver cancer data were also conducted by Kamo et al. (2011), Tondaet al. (2015) and so on.
We considered the Bayesian APC models described in Section 2 and conducted Bayesian inferencewith INLA as in Section 3. The random walk priors of order 1 or 2 (RW1 or RW2) were incorporated
for each of age, period and cohort eﬀects. The deviance information criterion (Spiegelhalter et al., 2002;Gelman et al., 2014)
DIC = −2 log p(y|ẑB , θ̂B ) + 2dDIC ,(4.1)
was computed for each combination of the priors, where ẑB and θ̂B are Bayes estimates of z andθ, respectively, and dDIC is the eﬀective number of parameters. Bayesian inference with INLA andcomputation of the DIC value were conducted by the R function inla of the package INLA.
The DIC values for some combinations of the priors are listed in Table 1. For the male data set, themodel with RW1 priors on age, period and cohort eﬀects gave the smallest DIC value, while for thefemale data set, the model with RW2 priors on age and period eﬀects and RW1 priors on cohort eﬀect
gave the smallest DIC value, The model with RW2 priors on age, period and cohort eﬀects fitted poorly
for both male and female data sets, so DIC values are not listed.
Figure 1 shows posterior means and 95% credible intervals of age, period and cohort eﬀects in themodel giving smallest DIC, fitted to male and female data. The age eﬀect was non-monotonic; the risk
was higher in infancy, getting lower until 15 years old, and getting higher with age. On the cohort
eﬀect, the risk was changed suddenly and became highest around the birth year of 1930. The suddenchange would explain the reason why the model with RW2 prior for the cohort eﬀect was poor, because
RW2 was too smooth. The period eﬀect was smaller than the age and cohort eﬀects.
It would be possible to provide prediction in future periods by extrapolation using the fitted
Bayesian APC models.
5Concluding remarks

We illustrated an application of the Bayesian age–period–cohort models to Japanese liver cancer mortality data. A sudden change of the cohort eﬀect was suggested by selecting a model with RW1 prior.
In our future works, we need to compare the performance of INLA with some other computation
methods MCMC and the empirical Bayes method with simple Laplace approximation in the Bayesian
APC models. Also, we should apply the described method to incidence and mortality data other thanliver cancer. Moreover, we intend to consider a Bayesian APC model with region eﬀects using cancer
data separated by prefecture.
Table 1: DIC values for some combinations of the priors
{αi }
RW1
RW2
RW2
RW2{β j }
RW1
RW1
RW2
RW1{γk }
RW1
RW1
RW1
RW2

Male
12050.15
12050.56
12052.4312054.91

Female8955.228
8957.157
8938.019
8953.035

0

20

40

60

80

age

1960

1980
period

2000

1920

1960

cohort2000

2040age60

80

0liver.female.cohort$mean−2−4−40

24f(cohort) (female)

420liver.female.period$mean

0liver.female.age$mean−2−41880−24f(period) (female)

24f(age) (female)

20liver.male.cohort$mean−4−220liver.male.period$mean−4−220−4−2liver.male.age$meanf(cohort) (male)

4f(period) (male)

4f(age) (male)

1960

1980

2000

1880

period

1920

1960

2000

cohortFigure 1: Posterior means and 95% credible intervals of age, period and cohort eﬀects.
Acknowledgment
This research was partly supported by JSPS KAKENHI Grant No. 26330042.
References[1] Carstensen, B. (2007). Age–period–cohort models for the Lexis diagram. Statistics in Medicine, 26, 3018–3045.
[2] Gelman, A., Hwang, J. and Vehtari, A. (2014). Understanding predictive information criteria for Bayesianmodels. Statistics and Computing, 24, 997–1016.
[3] Hills, M., Plummer, M. and Carstensen, B. (2014). A short introduction to R for epidemiology (Version 4).
http://bendixcarstensen.com/Epi/R-intro.pdf[4] Kamo, K., Tonda, T. and Satoh, K. (2011). Cancer mortality risk visualization on age-period plane. Proceedingsof the Institute of Statistical Mathematics, 59, 217–237 (in Japanese).
[5] Kupper, L. L., Janis, J. M., Karmous, A. and Greenberg, B. G. (1985). Statistical age–period–cohort analysis: areview and critique. Journal of Chronic Diseases, 38, 831–836.
[6] Martino, S. and Rue, H. (2010). Implementing approximate Bayesian inference using integrated nested Laplaceapproximation: a manual for the inla program. Department of Mathematical Sciences, NTNU, Norway.
[7] Nakamura, T. (1986). Bayesian cohort models for general cohort table analyses. Annals ofthe Institute ofStatistical Mathematics. 38, 353–370.
[8] Nielsen, B. (2015). apc: An R package for age–period–cohort analysis. The R Journal, 7, 52–64.
[9] Riebler, A., Held, L. and Rue, H. (2012). Estimation and extrapolation of time trends in registry data-borrowing
strength from related population. The Annals of Applied Statistics, 6, 304–333.
[10] Rue, H. and Held, L. (2005). Gaussian Markov Random Fields: Theory and Applications. Chapman and Hall.
[11] Rue, H., Martino, S. and Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian models byusing integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series
B, 71, 319–392.
[12] Schmid, V. J. and Held, L. (2007). Bayesian age–period–cohort modeling and prediction - BAMP. Journal ofStatistical Software, 21, 1–15.
[13] Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and van der Linde, A. (2002). Bayesian measures of model
complexity and fit (with discussion). Journal of the Royal Statistical Society, Series B, 64, 583–639.
[14] Tonda, T., Satoh, K. and Kamo, K. (2015). Detecting a local cohort eﬀect for cancer mortality data using avarying coeﬃcient model. Journal of Epidemiology, 25, 639–646.