###################################################
### chunk number 1: 
###################################################
library(Epi)


###################################################
### chunk number 2: 
###################################################
lung <- read.table( "../data/lung5-M.txt", header=T )
head(lung)
attach( lung )
table( A )
table( P )


###################################################
### chunk number 3: 
###################################################
round( tapply( Y, list(A,P), sum )/1000, 1 )


###################################################
### chunk number 4: 
###################################################
ap.1 <- glm( D ~ factor(A) + factor(P) + offset(log(Y)), 
             family=poisson, data=lung )
summary( ap.1 )


###################################################
### chunk number 5: 
###################################################
ap.0 <- glm( D ~ -1 +  factor(A) + factor(P) + offset(log(Y)), 
             family=poisson, data=lung )
summary( ap.0 )


###################################################
### chunk number 6: 
###################################################
ap.2 <- glm( D ~ factor(A) - 1 + relevel(factor(P),"1968") + offset(log(Y)), 
                 family=poisson, data=lung )


###################################################
### chunk number 7: 
###################################################
( ap.cf <- summary( ap.2 )$coef )


###################################################
### chunk number 8: 
###################################################
age.cf <- ap.cf[1:10,1:2]


###################################################
### chunk number 9: agerates
###################################################
par( mfrow=c(1,3) )
am <- seq(40,85,5)+2.5
plot( am, age.cf[,1] )
plot( am, exp(age.cf[,1]), log="y" )
plot( am, exp(age.cf[,1]), type="l", log="y" )


###################################################
### chunk number 10: agerates2
###################################################
matplot( am, cbind( exp(age.cf[,1]),
                    exp(age.cf[,1]-1.96*age.cf[,2]),
                    exp(age.cf[,1]+1.96*age.cf[,2]) ),
         type="l", log="y", lwd=c(3,1,1), lty=1, col="black" )

matplot( am, cbind( exp(age.cf[,1]),
                    exp(age.cf[,1]-1.96*age.cf[,2]),
                    exp(age.cf[,1]+1.96*age.cf[,2]) )*1000,
         xlab="Age", ylab="Rates of lung cancer in 1968 (per 1000 PY)",
         type="l", log="y", lwd=c(3,1,1), lty=1, col="black" )
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )

###################################################
### chunk number 11: 
###################################################
( RR.cf <- ap.cf[11:20,1:2] )


###################################################
### chunk number 12: 
###################################################
RR.cf <- rbind( RR.cf[1:5,], c(0,0), RR.cf[6:10,] )
RR.cf


###################################################
### chunk number 13: RR2
###################################################
matplot( as.numeric(levels(factor(P)))+2.5,
         cbind( exp(RR.cf[,1]),
                exp(RR.cf[,1]-1.96*RR.cf[,2]),
                exp(RR.cf[,1]+1.96*RR.cf[,2]) ),
         type="l", log="y", lwd=c(3,1,1), lty=1, col="black" )


###################################################
### chunk number 14: 
###################################################
ci.lin( ap.0, subset=c("A","1968") )


###################################################
### chunk number 15: 
###################################################
ci.lin( ap.0, subset=c("A","1968"), Exp=TRUE )


###################################################
### chunk number 16: 
###################################################
( cm.A <- cbind( diag( nlevels( factor(A) ) ), 1 ) )


###################################################
### chunk number 17: agerates3
###################################################
arates <- ci.lin( ap.0, subset=c("A","1968"), ctr.mat=cm.A, Exp=TRUE )[,5:7]
matplot( as.numeric( levels( factor(A) ) )+2.5, arates, 
         log="y", type="l", lwd=c(3,1,1), col="black", lty=1 )


###################################################
### chunk number 18: 
###################################################
cm.P <- rbind(0,diag( nlevels(factor(P))-1 ) )
cm.P


###################################################
### chunk number 19: 
###################################################
cm.Pref <- cm.P * 0
wh.col <- grep( "1968", levels(factor(P)) ) - 1 
cm.Pref[,wh.col] <- 1
cm.Pref


###################################################
### chunk number 20: RR3
###################################################
cm.P - cm.Pref
RR0 <- ci.lin( ap.0, subset="P", ctr.mat=cm.P-cm.Pref, Exp=TRUE )[,5:7]
matplot( as.numeric(levels(factor(P)))+2.5, RR0,
         type="l", log="y", lwd=c(3,1,1), lty=1, col="black" )


###################################################
### chunk number 21: 
###################################################
age.pt <- as.numeric(levels(factor(A)))+2.5  
 RR.pt <- as.numeric(levels(factor(P)))+2.5  
save( age.pt, arates,
       RR.pt, RR0, file="../data/age-per-est.Rdata" )


###################################################
### chunk number 22: 
###################################################
alim <- range( A ) + c(0,5) 
plim <- range( P ) + c(0,5)


###################################################
### chunk number 23: ageRR
###################################################
# Compute limits explicitly
rlim <- range(arates*10^5)*c(1/1.05,1.05)
RRlim <- 10^(log10(rlim)-ceiling(mean(log10(rlim))))         
# Determine relative width of plots
layout( rbind( c(1,2) ), widths=c(diff(alim),diff(plim)) )
# No space on the sides of the plots, only outer space
par( mar=c(4,0,1,0), oma=c(0,4,0,4), mgp=c(3,1,0)/1.5, las=1 )
matplot( as.numeric(levels(factor(A)))+2.5, arates*10^5,
         type="l", lwd=c(3,1,1), lty=1, col="black",
         log="y", xaxs="i", xlim=alim, xlab="Age", ylim=rlim )          
mtext( "Male lung cancer per 100,000", las=0, side=2, outer=T, line=2.5 )
matplot( as.numeric(levels(factor(P)))+2.5, RR0,
         type="l", lwd=c(3,1,1), lty=1, col="black",
         log="y", xlab="Period of follow-up", xlim=plim, yaxt="n", ylim=RRlim, ylab="" )
abline( h=1 )
points( 1968+2.5, 1, pch=1, lwd=3 )
axis( side=4 )
mtext( "Rate ratio", side=4, outer=T, las=0, line=2.5 )


#############################################
### Extras: Splines
#############################################
a.kn <- 5:7*10
a.bo <- c(40,85)
p.kn <- 1960
p.bo <- c(1950,1990)

########
### Illustration a spline basis
a.pr <- 40:90
matplot( a.pr, ns(a.pr, kn = a.kn, Bo = a.bo, int = T), type="l", lty=1, lwd=3)
matlines(a.pr, ns(a.pr, kn = a.kn, Bo = a.bo, int = F), type="l", lty=3, lwd=3)
abline(v=c(a.kn,a.bo))

### The model with splines
ap.s <- glm( D ~ ns(A,kn=a.kn,Bo=a.bo,int=T) +
                 ns(P,kn=p.kn,Bo=p.bo) - 1,
                 offset=log(Y/1000),
                 data=lung )
# The contrast