Skip to contents

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)

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:     230.12
## pD:            36.60
## DIC:          266.72
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.436        1.108        2.295
## period                          48.842      128.964      320.787
## cohort                          35.190       61.080      103.190
## 
## 
## 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)))
## [1] TRUE
print(model2)
## 
##  Model:
## age (rw2)  - period (rw2)  - cohort (rw2) model
## Deviance:     233.58
## pD:            35.72
## DIC:          269.30
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              1.090        3.024        6.847
## period                          16.052       41.275       91.027
## cohort                          24.825       46.878       83.398
## 
## 
## Markov Chains convergence checked succesfully using Gelman's R (potential scale reduction factor).
plot(model2)

model3<-bamp(cases, population, age="rw1", period=" ", cohort="rw2",
              periods_per_agegroup = 5)
## [1] TRUE
print(model3)
## 
##  Model:
## age (rw1) cohort (rw2) model
## Deviance:     275.82
## pD:            29.69
## DIC:          305.52
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.299        0.724        1.497
## cohort                          36.892       73.124      139.543
## 
## 
## 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:     230.43
## pD:            37.39
## DIC:          267.81
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.408        1.021        2.101
## period                          55.546      142.037      349.491
## cohort                          31.623       54.483       92.487
## 
## 
## 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:     229.64
## pD:            36.71
## DIC:          266.35
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.442        1.098        2.226
## period                          54.273      147.975      384.812
## cohort                          34.560       60.147       98.309
## 
## 
## Markov Chains convergence checked succesfully using Gelman's R (potential scale reduction factor).
plot(model5)