Simulating Age-Period-Cohort Data
Volker Schmid
2022-05-04
age=2*sqrt(seq(1,20,length=10))
age<- age-mean(age)
plot(age, type="l")
period=15:1
period[8:15]<-8:15
period<-period/5
period<-period-mean(period)
plot(period, type="l")
periods_per_agegroup=5
number_of_cohorts <- periods_per_agegroup*(10-1)+15
cohort<-rep(0,60)
cohort[1:15]<-(14:0)
cohort[16:30]<- (1:15)/2
cohort[31:60]<- 8
cohort<-cohort/10
cohort<-cohort-mean(cohort)
plot(cohort, type="l")
simdata<-apcSimulate(-10, age, period, cohort, periods_per_agegroup, 1e6)
print(simdata$cases)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 4 23 53 75 130 168 405 1121 3072
## [2,] 1 4 20 36 78 122 127 324 834 2253
## [3,] 0 3 11 36 67 86 127 210 634 1686
## [4,] 2 2 10 24 46 77 107 170 447 1294
## [5,] 1 4 7 25 38 66 81 103 333 898
## [6,] 1 1 9 16 34 54 74 112 240 679
## [7,] 0 1 6 11 28 53 82 98 173 499
## [8,] 0 3 8 10 20 55 46 89 134 381
## [9,] 0 1 5 18 26 44 75 113 156 426
## [10,] 0 1 7 23 38 72 97 138 178 424
## [11,] 2 4 6 19 53 82 113 164 216 548
## [12,] 1 11 16 19 52 111 163 193 307 538
## [13,] 0 2 17 41 71 126 200 272 317 613
## [14,] 1 8 16 40 81 177 260 348 446 666
## [15,] 1 7 26 43 98 197 327 495 607 750
simmod <- bamp(cases = simdata$cases, population = simdata$population, age = "rw1",
period = "rw1", cohort = "rw1", periods_per_agegroup =periods_per_agegroup)
##
## Automatic check procedure removed 1 Markov chain. Please check for convergence using checkConvergence() and maybe change your model settings (maybe add overdispersion).
## Warning: MCMC chains did not converge!
##
## WARNING! Markov Chains have apparently not converged! DO NOT TRUST THIS MODEL!
##
## Model:
## age (rw1) - period (rw1) - cohort (rw1) model
## Deviance: 142.57
## pD: 48.72
## DIC: 191.29
##
##
## Hyper parameters: 5% 50% 95%
## age 0.527 1.207 2.408
## period 13.246 26.223 46.296
## cohort 79.998 127.928 190.671
## Warning: MCMC chains did not converge!
## [1] FALSE
effects<-effects(simmod)
effects2<-effects(simmod, mean=TRUE)
#par(mfrow=c(3,1))
plot(age, type="l")
lines(effects$age, col="blue")
lines(effects2$age, col="green")
plot(period, type="l")
lines(effects$period, col="blue")
lines(effects2$period, col="green")
plot(cohort, type="l")
lines(effects$cohort, col="blue")
lines(effects2$cohort, col="green")
prediction<-predict_apc(simmod, periods=5, population=array(1e6,c(20,10)))
plot(prediction$cases_period[2,], ylim=range(prediction$cases_period),ylab="",pch=19)
points(prediction$cases_period[1,],pch="–",cex=2)
points(prediction$cases_period[3,],pch="–",cex=2)
for (i in 1:20)lines(rep(i,3),prediction$cases_period[,i])
plot(prediction$period[2,])
cov_p<-rnorm(15,period,.1)
simmod2 <- bamp(cases = simdata$cases, population = simdata$population, age = "rw1",
period = "rw1", cohort = "rw1", periods_per_agegroup =periods_per_agegroup,
period_covariate = cov_p)
##
## Model:
## age (rw1) - period (rw1) - cohort (rw1) model
## Deviance: 142.44
## pD: 48.65
## DIC: 191.09
##
##
## Hyper parameters: 5% 50% 95%
## age 0.517 1.197 2.407
## period 13.761 26.212 45.274
## cohort 78.635 127.149 191.685
##
##
## Markov Chains convergence checked succesfully using Gelman's R (potential scale reduction factor).
checkConvergence(simmod2)
## [1] TRUE
