# Clayton Hills, example 22.6 p221 

# [jh version 2008.02.24]

# DATA

cases=c(4,2,5,12,8,14);

pt=c(607.9,311.9,1272.1,878.1,888.9,667.5)

e=c(0,1,0,1,0,1)

# 2 age-category indicators *

a50=c(0,0,1,1,0,0)
a60=c(0,0,0,0,1,1)

# * could also use the gl function (below)
#   to set up variables that have a pattern
#   (or use the as.factor to create indicator -- 'dummy' --
#   variables for categorical variables). But you
#   have more control' if you set them up 'manually'

# no need to think of exposure as a `factor' since
# it is already 0/1 and can be included as a linear 
# (continuous) variable


#######################################
# to fit additive (ie ID Ratio) model
######################################

# make products with pt
# (you will need to fill in the blanks ... )

e.pt   =   e*pt 
a50.pt = a50*pt
a60.pt = ...

# for rate DIFFERENCE models, use POISSON variation
# and IDENTITY link. IDENTITY just means that you
# are asking that the (untransformed) 
# E[#cases] = linear predictor = [B0] + B1.X1 + B2.X2 + .. 
# you are not modelling some transform of E[#cases]
# (for the multiplicative model below, you will model
# a transsformation , ie the log, of E[#cases] as a
# function of the linear predictor 

glm(cases~-1 + pt+ e.pt + ... ,
    family=poisson(link="identity"))
    
fitted.cases = round(fitted(glm(cases~-1 + pt+ e.pt + ... ,
                     family=poisson(link="identity")) ),
                     1);

residuals = round(residuals(glm(cases~-1 + pt+ e.pt + ... ,
                     family=poisson(link="identity")) ),
                     1);

cbind(cases,fitted.cases,residuals)

# etc etc 

#############################################
# to fit multiplicative (ie ID Ratio) model
#############################################

# make levels for age and exposure
# (could also create manually, see e a50 a60 above)

e=gl(2,1,6,labels=c("1","0"));e
a=gl(3,2,6,labels=c("40","50","60"));a

# for rate RATIO (ie multiplicative) models, use 
# POISSON variation. By default, wehn you specify
# Poisson, the LOG link is used, i.e.,
# 
# log{E[#cases]} = linear predictor           + offset
# 
#                = B0 + B1.X1 + B2.X2 + ..    + offset  
# 
# This is equivalent to 
#
# E[#cases] = exp(linear predictor               + offset)
#
#           = exp(linear predictor)              * exp(offset)
#
#           = exp(B0)*exp(B1.X1)*exp(B2.X2)* ... * PT
#
#           = <--------- RATE------------------> * PT 

coef(glm(cases~e+a,family=poisson,offset=log(pt)))

exp(coef(glm(cases~e+a,family=poisson,offset=log(pt))))

exp(confint(glm(cases~e+a,family=poisson,offset=log(pt))))

fitted(glm(cases~e+a,family=poisson,offset=log(pt)))

residuals(glm(cases~e+a,family=poisson,offset=log(pt)))

cases

# etc etc