### R code from vignette source 'apc-tri-sol.rnw'

###################################################
### code chunk number 1: apc-tri-sol.rnw:19-23
###################################################
library( Epi )
ltri <- read.table( "../data/lung5-Mc.txt", header=T )
head( ltri )
ltri$S5 <- ltri$P5 - ltri$A5


###################################################
### code chunk number 2: Lexis-ill
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
Lexis.diagram( age=30+c(0,65), date=1938+c(0,65), coh.grid=TRUE )
with( ltri, text( Px, Ax, paste(D), cex=0.8, font=2 ) )
box()


###################################################
### code chunk number 3: apc-tri-sol.rnw:48-51
###################################################
ms <- glm( D ~ -1 + factor(A5) + factor(P5) + factor(S5) + offset(log(Y)),
           family=poisson, data=ltri )
summary( ms )$df


###################################################
### code chunk number 4: apc-tri-sol.rnw:56-59
###################################################
mc <- glm( D ~ -1 + factor(A5) + factor(P5) + factor(C5) + offset(log(Y)),
           family=poisson, data=ltri )
summary( mc )$df


###################################################
### code chunk number 5: parmest1
###################################################
par( mfrow=c(1,3) )
a.pt <- as.numeric( levels(factor(A5)) )
p.pt <- as.numeric( levels(factor(P5)) )
s.pt <- as.numeric( levels(factor(S5)) )
c.pt <- as.numeric( levels(factor(C5)) )

matplot( a.pt, ci.lin( ms, subset="A5", Exp=TRUE )[,5:7]/10^5,
         type="l", lty=1, lwd=c(3,1,1), col="black",
         xlab="Age", ylab="Rates", log="y" )
matlines( a.pt, ci.lin( mc, subset="A5", Exp=TRUE )[,5:7]/10^5,
          type="l", lty=1, lwd=c(3,1,1), col="blue" )

matplot( p.pt, rbind( c(1,1,1), ci.lin( ms, subset="P5",Exp=TRUE )[,5:7] ),
         type="l", lty=1, lwd=c(3,1,1), col="black",
         xlab="Period", ylab="RR", log="y" )
matlines( p.pt, rbind( c(1,1,1), ci.lin( mc, subset="P5",Exp=TRUE )[,5:7] ),
          type="l", lty=1, lwd=c(3,1,1), col="blue" )

matplot( s.pt, rbind(c(1,1,1),ci.lin( ms, subset="S5", Exp=TRUE )[,5:7]),
         type="l", lty=1, lwd=c(3,1,1), col="black",
         xlab="Cohort", ylab="RR", log="y" )
matlines( c.pt, rbind(c(1,1,1),ci.lin( mc, subset="C5", Exp=TRUE )[,5:7]),
          type="l", lty=1, lwd=c(3,1,1), col="blue" )


###################################################
### code chunk number 6: apc-tri-sol.rnw:112-115
###################################################
mt <- glm( D ~ -1 + factor(Ax) + factor(Px) + factor(Cx) + offset(log(Y)),
           family=poisson, data=ltri )
summary( mt )$df


###################################################
### code chunk number 7: parmest2
###################################################
par( mfrow=c(1,3) )
a.pt <- as.numeric( levels(factor(Ax)) )
p.pt <- as.numeric( levels(factor(Px)) )
c.pt <- as.numeric( levels(factor(Cx)) )

matplot( a.pt, ci.lin( mt, subset="Ax", Exp=TRUE )[,5:7]/10^5,
         type="l", lty=1, lwd=c(3,1,1), col="black",
         xlab="Age", ylab="Rates", log="y" )

matplot( p.pt, rbind( c(1,1,1), ci.lin( mt, subset="Px",Exp=TRUE )[,5:7] ),
         type="l", lty=1, lwd=c(3,1,1), col="black",
         xlab="Period", ylab="RR", log="y" )

matplot( c.pt, rbind(c(1,1,1),ci.lin( mt, subset="Cx", Exp=TRUE )[,5:7]),
         type="l", lty=1, lwd=c(3,1,1), col="black",
         xlab="Cohort", ylab="RR", log="y" )


###################################################
### code chunk number 8: apc-tri-sol.rnw:148-149
###################################################
summary( mt )$deviance


###################################################
### code chunk number 9: apc-tri-sol.rnw:155-156
###################################################
table( up, P5-A5-C5 )


###################################################
### code chunk number 10: apc-tri-sol.rnw:164-172
###################################################
m.up <- glm( D ~ -1 + factor(A5) + factor(P5) + factor(S5) + offset(log(Y)),
             family=poisson, data=subset(ltri,up==1) )
summary( m.up )$deviance
m.lo <- glm( D ~ -1 + factor(A5) + factor(P5) + factor(S5) + offset(log(Y)),
             family=poisson, data=subset(ltri,up==0) )
summary( m.lo )$deviance
summary( m.lo )$deviance + summary( m.up )$deviance
summary( mt )$deviance


###################################################
### code chunk number 11: parmest3
###################################################
par( mfrow=c(1,3) )
a.pt <- as.numeric( levels(factor(Ax)) )
p.pt <- as.numeric( levels(factor(Px)) )
c.pt <- as.numeric( levels(factor(Cx)) )

a5.pt <- as.numeric( levels(factor(A5)) )
p5.pt <- as.numeric( levels(factor(P5)) )
s5.pt <- as.numeric( levels(factor(S5)) )

matplot( a.pt, ci.lin( mt, subset="Ax", Exp=TRUE )[,5:7]/10^5,
         type="l", lty=1, lwd=c(2,1,1), col=gray(0.7),
         xlab="Age", ylab="Rates", log="y" )
matpoints( a5.pt, ci.lin( m.up, subset="A5", Exp=TRUE )[,5:7]/10^5,
           pch=c(16,3,3), col="blue" )
matpoints( a5.pt, ci.lin( m.lo, subset="A5", Exp=TRUE )[,5:7]/10^5,
           pch=c(16,3,3), col="red" )

matplot( p.pt, rbind( c(1,1,1), ci.lin( mt, subset="Px",Exp=TRUE )[,5:7] ),
         type="l", lty=1, lwd=c(2,1,1), col=gray(0.7),
         xlab="Period", ylab="RR", log="y" )
matpoints( p5.pt[-1], ci.lin( m.up, subset="P5", Exp=TRUE )[,5:7],
           pch=c(16,3,3), col="blue" )
matpoints( p5.pt[-1], ci.lin( m.lo, subset="P5", Exp=TRUE )[,5:7],
           pch=c(16,3,3), col="red" )

matplot( c.pt, rbind(c(1,1,1),ci.lin( mt, subset="Cx", Exp=TRUE )[,5:7]),
         type="l", lty=1, lwd=c(2,1,1), col=gray(0.7),
         xlab="Cohort", ylab="RR", log="y" )
matpoints( s5.pt[-1], ci.lin( m.up, subset="S5", Exp=TRUE )[,5:7],
           pch=c(16,3,3), col="blue" )
matpoints( s5.pt[-1], ci.lin( m.lo, subset="S5", Exp=TRUE )[,5:7],
           pch=c(16,3,3), col="red" )


###################################################
### code chunk number 12: apc-tri-sol.rnw:236-244
###################################################
library( splines )
mspl <- glm( D ~ -1 + ns(Ax,df=7,intercept=T)
                    + ns(Px,df=6,intercept=F)
                    + ns(Cx,df=6,intercept=F) + offset(log(Y)),
             family=poisson, data=ltri )
summary( mspl )
summary( mt )$deviance - summary( mspl )$deviance
summary( mt )$df       - summary( mspl )$df


###################################################
### code chunk number 13: parmest4
###################################################
pspl <- predict( mspl, type="terms", se.fit=TRUE )
str(pspl)
a.ord <- order( ltri$Ax )
p.ord <- order( ltri$Px )
c.ord <- order( ltri$Cx )

par( mfrow=c(1,3) )
matplot( Ax[a.ord], exp(cbind( pspl$fit[,1], pspl$se.fit[,1] )[a.ord,] %*% ci.mat())*10^5,
         type="l", lty=1, lwd=c(2,1,1), col=gray(0.2),
         xlab="Age", ylab="Rates", log="y" )

matplot( Px[p.ord], exp(cbind( pspl$fit[,2], pspl$se.fit[,2] )[p.ord,] %*% ci.mat()),
         type="l", lty=1, lwd=c(2,1,1), col=gray(0.2),
         xlab="Period", ylab="RR", log="y" )

matplot( Cx[c.ord], exp(cbind( pspl$fit[,3], pspl$se.fit[,3] )[c.ord,] %*% ci.mat()),
         type="l", lty=1, lwd=c(2,1,1), col=gray(0.2),
         xlab="Cohort", ylab="RR", log="y" )