Skip to contents

Bayesian Age-Period-Cohort Modeling

Volker Schmid

2026-07-06

Source: vignettes/largevignettes/modeling.Rmd
modeling.Rmd

Data example

BAMP includes a data example.

data(apc)
plot(cases[,1],type="l",ylim=range(cases), ylab="cases", xlab="year", main="cases per age group")
for (i in 2:8)lines(cases[,i], col=i)

APC model with random walk first order prior 

model1 <- bamp(cases, population, age="rw1", period="rw1", cohort="rw1",
              periods_per_agegroup = 5)
## 
## Automatic check procedure removed 1 Markov chain. Please check for convergence using checkConvergence() and maybe change your model settings (maybe add overdispersion).

bamp() automatically performs a check for MCMC convergence using Gelman and Rubin’s convergence diagnostic. We can manually check the convergence again:

## [1] TRUE

Now we have a look at the model results. This includes estimates of smoothing parameters and deviance and DIC:

print(model1)
## 
##  Model:
## age (rw1)  - period (rw1)  - cohort (rw1) model
## Deviance:     231.16
## pD:            36.79
## DIC:          267.96
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.360        0.930        1.970
## period                          67.829      195.343      602.609
## cohort                          34.501       59.455       98.374
## 
## 
## Markov Chains convergence checked succesfully using Gelman's R (potential scale reduction factor).

We can plot the main APC effects using point-wise quantiles:

plot(model1)

More quantiles are possible:

plot(model1, quantiles = c(0.025,0.1,0.5,0.9,0.975))

APC model with random walk second order prior 

model2 <- bamp(cases, population, age="rw2", period="rw2", cohort="rw2",
              periods_per_agegroup = 5,
              mcmc.options=list("number_of_iterations"=200000, "burn_in"=100000, "step"=50, "tuning"=500),
              hyperpar=list("age"=c(1,.5), "period"=c(1,0.05), "cohort"=c(1,0.05)))
## Warning: MCMC chains did not converge!
## Warning: MCMC chains did not converge!
## [1] FALSE
print(model2)
## 
## WARNING! Markov Chains have apparently not converged! DO NOT TRUST THIS MODEL!
## 
##  Model:
## age (rw2)  - period (rw2)  - cohort (rw2) model
## Deviance:     233.99
## pD:            36.84
## DIC:          270.84
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              1.076        2.921        6.751
## period                          16.195       41.631       90.001
## cohort                          23.438       44.912       82.160
plot(model2)

model3<-bamp(cases, population, age="rw1", period=" ", cohort="rw2",
              periods_per_agegroup = 5)
## 
## Automatic check procedure removed 1 Markov chain. Please check for convergence using checkConvergence() and maybe change your model settings (maybe add overdispersion).
## [1] TRUE
print(model3)
## 
##  Model:
## age (rw1) cohort (rw2) model
## Deviance:     276.39
## pD:            30.08
## DIC:          306.47
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.280        0.722        1.500
## cohort                          38.858       73.503      136.295
## 
## 
## Markov Chains convergence checked succesfully using Gelman's R (potential scale reduction factor).
plot(model3)

(model4<-bamp(cases, population, age="rw1", period="rw1", cohort="rw1",
             cohort_covariate = cov_c, periods_per_agegroup = 5))
## 
##  Model:
## age (rw1)  - period (rw1)  - cohort (rw1) model
## Deviance:     231.26
## pD:            36.92
## DIC:          268.17
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.351        0.906        1.896
## period                          67.027      202.564      637.852
## cohort                          34.056       58.626       98.566
## 
## 
## Markov Chains convergence checked succesfully using Gelman's R (potential scale reduction factor).
plot(model4)

(model5<-bamp(cases, population, age="rw1", period="rw1", cohort="rw1",
             period_covariate = cov_p, periods_per_agegroup = 5))
## 
##  Model:
## age (rw1)  - period (rw1)  - cohort (rw1) model
## Deviance:     231.17
## pD:            36.91
## DIC:          268.08
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.369        0.913        1.938
## period                          66.587      199.151      600.304
## cohort                          34.579       59.617       97.261
## 
## 
## Markov Chains convergence checked succesfully using Gelman's R (potential scale reduction factor).
plot(model5)