#====================================================
# Illustration of Poisson regression for grouped data
# Example: Lung cancer in Danish cities
#====================================================

# Write the data into vectors:
city=rep(1:4,each=5)
age=rep(1:5,times=4)
cases=c(11,11,11,10,11,13,6,15,10,12,4,8,7,11,9,5,7,10,14,8)
number=c(3059,800,710,581,509,2879,1083,923,834,634,3142,1050,895,702,535,2520,878,839,631,539)

# Create data frame
lungcancer=data.frame(city,age,number,cases)
lungcancer

# Fit Poisson GLM model with log link and city and age group 
# as factors and no intercept:
fit1=glm(cases~offset(log(number))+factor(age)+factor(city)-1, family=poisson,data=lungcancer)
summary(fit1)

# Plot of baseline rate (i.e. effect of age for city=1):
age.m=c(47.5,57.5,62.5,67.5,72.5)   # mean age in each of the five age intervals
base.rate=exp(fit1$coef[1:5])
plot(age.m,base.rate*1000,xlab="Age",ylab="Rate per 1000",
     xlim=c(45,75),ylim=c(0,25),type='b',col='red')

# Print rate ratios of cities 2, 3 and 4 (relative to city 1, for fixed age group)
exp(fit1$coef[6:8])


# Test for effect of city
fit2=glm(cases~offset(log(number))+factor(age),family=poisson)
summary(fit2)
anova(fit2,fit1,test="LRT")