ReproducingMMNN2015.R
In apc: Age-Period-Cohort Analysis

### R code from vignette source 'ReproducingMMNN2015.Rnw'

###################################################
### code chunk number 1: ReproducingMMNN2015.Rnw:94-96
###################################################
library(apc)
data <- data.asbestos()


###################################################
### code chunk number 2: ReproducingMMNN2015.Rnw:107-108
###################################################
apc.plot.data.all(data)


###################################################
### code chunk number 3: ReproducingMMNN2015.Rnw:119-120
###################################################
apc.plot.data.sums(data)


###################################################
### code chunk number 4: ReproducingMMNN2015.Rnw:123-137
###################################################
data.sums <- apc.data.sums(data)
par(mfrow=c(2,2),oma=c(0,0,2,0),mar=c(4,4,2,0)+0.1)
plot(seq(25,89),data.sums$sums.age,
    main="(a) sums by age",
    type="l",xlab="age",ylab="age data sum")
plot(seq(1967,2007),data.sums$sums.per,
    main="(b) sums by period",
    type="l",xlab="period",ylab="period data sum")
plot(seq(1878,1982),data.sums$sums.coh,
    main="(c) sums by cohort",
    type="l",xlab="cohort",ylab="cohort data sum")
apc.plot.data.within(data,plot.type="pwc",
    thin=5,type="l",main="(d)",lty=1,legend=FALSE)
title("Figure 1",outer=TRUE)    


###################################################
### code chunk number 5: ReproducingMMNN2015.Rnw:144-145
###################################################
apc.fit.table(data,"poisson.response")[1:4,1:6]


###################################################
### code chunk number 6: ReproducingMMNN2015.Rnw:153-154
###################################################
fit.apc <- apc.fit.model(data,"poisson.response","APC")


###################################################
### code chunk number 7: ReproducingMMNN2015.Rnw:158-159
###################################################
apc.plot.fit.all(fit.apc)


###################################################
### code chunk number 8: ReproducingMMNN2015.Rnw:171-172
###################################################
apc.plot.fit.pt(fit.apc,main="Figure 3 - approximately")


###################################################
### code chunk number 9: ReproducingMMNN2015.Rnw:180-181
###################################################
fit.ac <- apc.fit.model(data,"poisson.response","AC")


###################################################
### code chunk number 10: ReproducingMMNN2015.Rnw:203-206
###################################################
forecast.16 <- apc.forecast.ac(fit.ac,sum.per.by.coh=c(42,89))
forecast.30 <- apc.forecast.ac(fit.ac,sum.per.by.coh=c(42,75))
forecast.45 <- apc.forecast.ac(fit.ac,sum.per.by.coh=c(42,60))


###################################################
### code chunk number 11: ReproducingMMNN2015.Rnw:213-214
###################################################
data.sum.per <- apc.data.sums(data.asbestos())$sums.per


###################################################
### code chunk number 12: ReproducingMMNN2015.Rnw:231-242
###################################################
plot(seq(1967,2007),data.sum.per,xlim=c(1967,2047),ylim=c(0,3500),
    xlab="period",ylab="number of cases",
    main="Figure 6")
apc.polygon(forecast.16$response.forecast.per.by.coh,
    2007,TRUE,TRUE,col.line=1)
apc.polygon(forecast.30$response.forecast.per.by.coh,
    2007,TRUE,TRUE,col.line=2)
apc.polygon(forecast.45$response.forecast.per.by.coh,
    2007,TRUE,TRUE,col.line=3)
legend("topleft",legend=c("1966","1952","1937"),
    col=c(1,2,3),lty=1)


###################################################
### code chunk number 13: ReproducingMMNN2015.Rnw:255-258
###################################################
data.1991 <- apc.data.list.subset(data.asbestos(),0,0,0,16,0,0)
fit.ac.1991 <- apc.fit.model(data.1991,"poisson.response","AC")
forecast.1991 <- apc.forecast.ac(fit.ac.1991)


###################################################
### code chunk number 14: ReproducingMMNN2015.Rnw:269-277
###################################################
data.2001 <- apc.data.list.subset(data.asbestos(),0,0,0,6,0,0)
fit.ac.2001 <- apc.fit.model(data.2001,"poisson.response","AC")
forecast.2001 <- apc.forecast.ac(fit.ac.2001)
data.2006 <- apc.data.list.subset(data.asbestos(),0,0,0,1,0,0)
fit.ac.2006 <- apc.fit.model(data.2006,"poisson.response","AC")
forecast.2006 <- apc.forecast.ac(fit.ac.2006)
fit.ac.2007 <- apc.fit.model(data.asbestos(),"poisson.response","AC")
forecast.2007 <- apc.forecast.ac(fit.ac.2007)


###################################################
### code chunk number 15: ReproducingMMNN2015.Rnw:282-297
###################################################
plot(seq(1967,2007),data.sum.per,xlim=c(1967,2047),ylim=c(0,3500),
    xlab="period",ylab="number of cases",
    main="Figure 7")
apc.polygon(forecast.2007$response.forecast.per.ic,
    2007,TRUE,TRUE,col.line=1)
apc.polygon(forecast.2007$response.forecast.per   ,
    2007,FALSE    ,col.line=2)
apc.polygon(forecast.2006$response.forecast.per   ,
    2006,FALSE    ,col.line=3)
apc.polygon(forecast.2001$response.forecast.per   ,
    2001,FALSE    ,col.line=4)
apc.polygon(forecast.1991$response.forecast.per   ,
    1991,FALSE    ,col.line=5)
legend("topleft",legend=c("2007ic","2007","2006","2001","1991"),
    col=c(1,2,3,4,5),lty=1)


###################################################
### code chunk number 16: ReproducingMMNN2015.Rnw:305-307
###################################################
forecast.2007$response.forecast.per[10:14,1]
forecast.2007$response.forecast.per.ic[10:14,1]


###################################################
### code chunk number 17: ReproducingMMNN2015.Rnw:323-329
###################################################
data.coh.1966 <- apc.data.list.subset(data.asbestos(),0,0,0,0,0,16)
fit.ac.coh.1966 <- apc.fit.model(data.coh.1966,"poisson.response","AC")
forecast.coh.1966 <- apc.forecast.ac(fit.ac.coh.1966)
data.age.35 <- apc.data.list.subset(data.asbestos(),10,0,0,0,0,0)
fit.ac.age.35 <- apc.fit.model(data.age.35,"poisson.response","AC")
forecast.age.35 <- apc.forecast.ac(fit.ac.age.35,sum.per.by.coh=c(42,89))


###################################################
### code chunk number 18: ReproducingMMNN2015.Rnw:339-341
###################################################
fit.apc.1966 <- apc.fit.model(data.coh.1966,"poisson.response","APC")
forecast.apc.1966   <- apc.forecast.apc(fit.apc.1966)


###################################################
### code chunk number 19: ReproducingMMNN2015.Rnw:351-370
###################################################
plot(seq(1967,2007),data.sum.per,xlim=c(1967,2047),ylim=c(0,3500),
    xlab="period",ylab="number of cases",
    main="Figure 8")
apc.polygon(forecast.16$response.forecast.per.by.coh    ,
    2007,FALSE,col.line=1)
apc.polygon(forecast.coh.1966$response.forecast.per     ,
    2007,FALSE,col.line=2)
apc.polygon(forecast.age.35$response.forecast.per.by.coh,
    2007,FALSE,col.line=3)
apc.polygon(forecast.16$response.forecast.per.by.coh.ic ,
    2007,FALSE,col.line=4)
apc.polygon(forecast.apc.1966$response.forecast.per     ,
    2007,FALSE,col.line=5)
legend("topleft",legend=c("est: full, fore: coh<=1966",
                          "est: coh<=1966, fore: coh<=1966",
                          "est: age>=35, fore: coh<=1966",
                          "est: full, fore: coh<=1966, ic",
                          "est: apc, I0: coh<=1966"),
    col=c(1,2,3,4,5),lty=1)


###################################################
### code chunk number 20: ReproducingMMNN2015.Rnw:377-381
###################################################
plot(seq(1967,2007),data.sum.per,xlim=c(1967,2047),ylim=c(0,3500),
    xlab="period",ylab="number of cases",
    main="Figure 8")
apc.polygon(forecast.16$response.forecast.per.by.coh.ic,2007,TRUE,TRUE)