library( Epi )
R <- c(6.1,  24,  72, 147, 192, 1.1, 11, 49, 108, 212)
D <- c( 32, 104, 206, 186, 102, 2  , 12, 28,  28,  31)
Y <- D / R # risk time in units of 10^4 PY
smk <- factor(rep(1:2, each = 5), labels = c("Smoke", "no-Sm") )
age <- factor(rep(seq(35, 75, 10), 2))
data.frame(D, Y, R, D/Y, age, smk)

# model with log link
ma <- glm(cbind(D, Y) ~ age + smk, family = poisreg)
mi <- update(ma, . ~ . + age:smk) # add the multiplicative interaction
anova(ma, mi, test = "Chisq")

# with identity link
aa <- glm(cbind(D, Y) ~ age + smk, family = poisreg(link = identity))
ai <- update(aa, . ~ . + age:smk ) # add the additive interaction
anova(aa, ai, test = "Chisq")

# multiplicative interaction
mI <- update(mi, . ~ -1 + age / smk)
ci.exp(mI, subset = "smk")
round(1 / ci.exp(mI, subset = "smk"), 2)

# additive interaction
aI <- update(ai, . ~ -1 + age / smk)
ci.exp(aI, Exp = FALSE)
round(- ci.exp(aI, Exp = FALSE, subset = "smk"), 2)

# forest plots
par(mfrow = c(1,2))
plotEst(1 / ci.exp(mI, subset = "smk"),
        xlog = TRUE, xlab = "smoking RR")
abline(v = c(1, exp(-coef(ma)[6])))
plotEst(-ci.exp(aI, Exp = FALSE, subset = "smk"),
        xlab = "smoking RD")
abline(v = c(0, -coef(aa)[6]))

# forest plots again, niceified
par(mfrow = c(1,2))
plotEst(1 / ci.exp(mI, subset = "smk"),
        xlog = TRUE, xlab = "smoking RR", xlim = c(0.5, 20))
abline(v = c(1, exp(-coef(ma)[6])), col = 1:2)
plotEst(-ci.exp(aI, Exp = FALSE, subset = "smk"),
        xlab = "smoking RD", xlim = c(-50, 100))
abline(v = c(0, -coef(aa)[6]), col = 1:2)
