# Congdon P., pag. 36, Example 2.1  Area mortality comparisons
#
# Short description.
# Application of the Poisson to compare mortality between areas after 
# standardizing for age.
# Silcocks (1994) presents data on male myeloid leukaemia deaths (1989) in 
# Derby, denoted r_1, and in the remainder of the Trent region of
# England, denoted r_2.
# The goal is comparing mortality between Derby, the index area, and
# the entire Trent region defined as the standard population. 
#
# Suppose r[i,j] and n[i,j] (i=1,...,6, j=1,2) denote the observed deaths
# and, respectively, the populations in area j for age groups i.
# If death rates in the standard population are e[i], the expected 
# deaths E in the index population are just sum_i(e[i]*n[i,1]).
# If actual deaths are equal, or nearly so, to expected deaths then 
# the mortality in the index area appears comparable to that in the 
# standard population.
#
# Model A)
# 
# A frequently used model assumes index area deaths r_1 are Poisson
# with mean E*rho where rho=1 if the standard and the index 
# death rates are the same (rho*100 = SMR, a standard mortality
# ratio, equalling 100 or nearly so when mortality is similar in 
# standard and index populations).
#
# Model B)
#
# While the deaths rates in the standard population are usually 
# assumed fixed, sometime they may be considered subject to sampling
# variation.
# Then, r[i,j] are considered Poisson with means lambda[i,j]=
# theta[i,j]*n[i,j].
# The expected deaths E in the index area are then calculated by
# E=sum_i(pstar[i](lambda[i,1]+lambda[i,2])),
# where pstar[i]=n[i,1]/(n[i,1]+n[i,2]) is the share of the total 
# standard population in age group i located in the index area.
# If the index population is a relatively large share of the standard 
# population there will be covariance between r_1=sum(r[,1]) and E
# and the credible interval of the SMR will be narrower than if the 
# expected deaths are treated as fixed.
# In program 2.15 this approach is denoted `standardization method 2'. 


# The program for A) is a simplification of the Program 2.15 
# Trent Leukaemia Mortality which fits model B).
# Program 2.15 can be downloaded from 
# http://www.mrc-bsu.cam.ac.uk/bugs/weblinks/webresource.shtml
# 
Program 2.15
Area Mortality Model 
 r[6,2],n[6,2],# deaths r, populations n; index popn. (1),
                     # standard popn. excluding index popn (2), by age (6 age groups)
 D,E,e[6],   # D=sum(r[,1]), total index deaths, E = sum(e), expected index deaths
 pstar[6],   # propn. of standard popn in age group i living in index district
 lambda[6,2],# Poisson deaths and underlying rates
 q.c[6],q,      # population shares weighted by deaths

{for (i in 1:Nage) {
# standardisation method 1 (Expected based on whole Trent region)
  pstar[i] <- n[i,1]/(n[i,1]+n[i,2]);e[i] <- pstar[i]*(lambda[i,1]+lambda[i,2]);
# standardisation method 2 (Expected based on Trent region excluding Derby)
#  pstar[i] <- n[i,1]/n[i,2];e[i] <- pstar[i]*lambda[i,2];
# weighted pstar
   q.c[i] <- pstar[i]*lambda[i,1];
# stochastic variability in both standard and index popns
                      for (j in 1:Narea) { r[i,j] ~ dpois(lambda[i,j]);
# underlying death rates
                               theta[i,j] ~ dgamma(0.001,0.001)
                               lambda[i,j] <- theta[i,j]*n[i,j];}}
# stochastic deaths and expected death totals
  D <- sum(lambda[,1]);
  E <- sum(e[]);
# SMR
  R <- D/E;
# weighted fraction of standard population represented by index popn
  q <- sum(q.c[])/D;
# covariance between index and expected deaths
  cov.D.E <- q*D }

Data
list(r=structure(.Data=c(1,0,1,0,3,7,9,23,7,16,9,30),.Dim=c(6,2)),
n=structure(.Data=c(29.9,125.6,28.1,264.8,206.5,514.0,103.2,405.7,
38.1,150.3,22.3,86.7),.Dim=c(6,2)),
Nage=6,Narea=2)


Model A)

#r[6,2], n[6,2],   # deaths r, populations n;
                   # index popn. (1), standard popn. including index popn (2),
                   # by age (6 age groups)

#D,E,e[6],   # D=sum(r[,1]), total index deaths, E = sum(e), expected index deaths

{
D~dpois(lambda)
rho~dgamma(0.001,0.001)
lambda<-E*rho
}


# Data for A)
# From the data supplied for model B:
# list( r=structure(.Data=c(1,0,1,0,3,7,9,23,7,16,9,30),.Dim=c(6,2)),
# n=structure(.Data=c(29.9,125.6,28.1,264.8,206.5,514.0,103.2,405.7,38.1,
# 150.3,22.3,86.7),.Dim=c(6,2))),
# D and E are calculated as 
# D<-sum(r[,1])
# for(i in 1:6) e[i] <- r[i,]/n[i,] 
# E<-sum(e[])
# so that the data for A) are the following:

list(E=22.38,D=30)