library(Epi)

# Male bladder cancer in Italy
data(blcaIT)
str(blcaIT)
bl <- transform( blcaIT, A=age+2.5, P=period+2.5 )
head( bl )

# A simple drift model
lAP <- glm( D ~ -1 + factor(A) + P,
            offset=log(Y/10^5),
            family=poisson,
            data=bl )
round( ci.lin( lAP ), 3 )

# Model with quadratic effect of period
qAP <- glm( D ~ -1 + factor(A) + P + I(P^2),
            offset=log(Y/10^5),
            family=poisson,
            data=bl )

# Prediction points for age
A.pr <- unique(bl$A)
# Prediciton points for period
P.pr <- 1955:1980
# Period reference - ALWAYS in a variable, you may want to change it later
P.ref <- 1970
# Contrast matrix (design matrix for age effects
CA <- model.matrix( ~ -1 + factor(A.pr) )
# Contrast matrices for period effects
CP     <- cbind( P.pr , P.pr^2 )
CP.ref <- cbind( P.ref, P.ref^2 )

# The rates in 1970 --- note how the CP.ref is expandend to match the CA.
R1970 <- ci.lin( qAP, subset=c("A","P"), ctr.mat=cbind(CA,CP.ref[rep(1,nrow(CA)),]), Exp=TRUE )[,5:7]
# The RR relative to 1970
RR.P  <- ci.lin( qAP, subset=      "P" , ctr.mat=      CP-CP.ref[rep(1,nrow(CP)),] , Exp=TRUE )[,5:7]

# Now plot the rates and the RRs
par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )
matplot( A.pr, R1970, type="l", col="black", lty=1, log="y", lwd=c(2,1,1) )
matplot( P.pr, RR.P , type="l", col="black", lty=1, log="y", lwd=c(2,1,1) )
abline(h=1,v=1970,col="gray")