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:     231.25
## pD:            36.95
## DIC:          268.20
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.339        0.891        1.959
## period                          67.892      195.109      589.123
## cohort                          34.303       59.298       99.935
## 
## 
## 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!
## [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:     234.90
## pD:            37.60
## DIC:          272.50
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              1.059        2.956        6.431
## period                          16.537       41.776       91.522
## cohort                          22.936       43.661       79.614
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:     276.60
## pD:            30.10
## DIC:          306.70
## 
## 
##  Hyper parameters:                 5%           50%          95%         
## age                              0.286        0.742        1.510
## cohort                          38.794       74.862      142.919
## 
## 
## 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))
plot(model4)

(model5<-bamp(cases, population, age="rw1", period="rw1", cohort="rw1",
             period_covariate = cov_p, periods_per_agegroup = 5))
plot(model5)