#### Ice: non-parametric smoothing in an age-cohort model

Breslow and Clayton (1993) analyse breast cancer rates in Iceland by year of birth (K = 11 cohorts from 1840-1849 to 1940-1949) and by age (J =13 groups from 20-24 to 80-84 years). Due to the number of empty cells we consider a single indexing over I = 77 observed number of cases, giving data of the following form.

In order to pull in the extreme risks associated with small birth cohorts, Breslow and
Clayton first consider the exchangeable model

casesi   ~   Poisson(μi)
log μi   =   log person-yearsi + αagei + βyeari
βk   ~   Normal( 0, τ )

Autoregressive smoothing of relative risks

They then consider the alternative approach of smoothing the rates for the cohorts by assuming an auto-regressive model on the β's, assuming the second differences are independent normal variates. This is equivalent to a model and prior distribution
casesi   ~   Poisson(μi)
log μi   =   log person-yearsi + αagei + βyeari
β1   ~   Normal( 0, 0.000001τ )
β2 | β1   ~   Normal( 0, 0.000001τ )
βk | β1,...,k-1   ~   Normal( 2 βk-1- βk-2, τ ) k > 2

We note that β1 and β2 are given "non-informative" priors, but retain a τ term in order to provide the appropriate likelihood for τ.

For computational reasons Breslow and Clayton impose constraints on their random effects βk in order that their mean and linear trend are zero, and counter these constraints by introducing a linear term b x yeari and allowing unrestrained estimation of αj. Since we allow free movement of the β's we dispense with the linear term, and impose a "corner" constraint α1 =0 .
model
{
for (i in 1:I) {
cases[i] ~ dpois(mu[i])
log(mu[i]) <- log(pyr[i]) + alpha[age[i]] + beta[year[i]]
}
betamean[1] <- 2 * beta[2] - beta[3]
Nneighs[1] <- 1
betamean[2] <- (2 * beta[1] + 4 * beta[3] - beta[4]) / 5
Nneighs[2] <- 5
for (k in 3 : K - 2) {
betamean[k] <- (4 * beta[k - 1] + 4 * beta[k + 1]- beta[k - 2] - beta[k + 2]) / 6
Nneighs[k] <- 6
}
betamean[K - 1] <- (2 * beta[K] + 4 * beta[K - 2] - beta[K - 3]) / 5
Nneighs[K - 1] <- 5
betamean[K] <- 2 * beta[K - 1] - beta[K - 2]
Nneighs[K] <- 1
for (k in 1 : K) {
betaprec[k] <- Nneighs[k] * tau
}
for (k in 1 : K) {
beta[k] ~ dnorm(betamean[k], betaprec[k])
logRR[k] <- beta[k] - beta[5]
tau.like[k] <- Nneighs[k] * beta[k] * (beta[k] - betamean[k])
}
alpha[1] <- 0.0
for (j in 2 : Nage) {
alpha[j] ~ dnorm(0, 1.0E-6)
}
d <- 0.0001 + sum(tau.like[]) / 2
r <- 0.0001 + K / 2
tau ~ dgamma(r, d)
sigma <- 1 / sqrt(tau)
}

##### Data
``` list(age = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6,             6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11,             11, 11, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13),   year = c(6, 7, 8, 9, 10, 11, 6, 7, 8, 9, 10, 11, 5, 6, 7, 8, 9, 10, 5, 6, 7, 8, 9, 10, 4, 5, 6, 7, 8, 9, 4, 5, 6, 7,             8, 9, 3, 4, 5, 6, 7, 8, 3, 4, 5, 6, 7, 8, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4,             5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5),   cases = c(2, 0, 1, 1, 1, 2, 0, 2, 1, 1, 5, 5, 1, 1, 3, 7, 12, 10, 6, 11, 9, 14, 20, 14, 7, 14, 22, 25, 29, 37, 21, 11, 29, 33,             57, 24, 15, 8, 22, 27, 38, 52, 10, 15, 22, 26, 47, 31, 8, 11, 17, 23, 31, 38, 8, 10, 24, 30, 53, 26, 5, 3, 10, 18,             22, 30, 1, 7, 11, 26, 32, 17, 5, 8, 17, 32, 31),   pyr = c(41380, 43650, 49810, 58105, 57105, 76380, 39615, 42205, 48315, 56785, 55965, 33955, 29150, 38460,            40810, 47490, 55720, 55145, 27950, 37375, 39935, 46895, 54980, 27810, 25055, 27040, 36400, 39355,         46280, 54350, 24040, 26290, 35480, 38725, 45595, 25710, 22890, 23095, 25410, 34420, 37725, 44740,          21415, 21870, 24240, 33175, 36345, 21320, 17450, 19765, 20255, 22760, 31695, 34705, 15350, 17720,          18280, 20850, 29600, 15635, 9965, 12850, 15015, 15725, 18345, 26400, 8175, 11020, 13095, 14050,         16480, 10885, 7425, 10810, 12260, 14780, 13600),I = 77, Nage = 13, K = 11) ```
##### Inits for chain 1
``` list(tau=1,alpha=c(NA,0,0,0,0,0,0,0,0,0,0,0,0),beta =c(0.05,0.1,0,0,0,0,0,0,0,0,0)) ```

##### Inits for chain 2
``` list(tau=1,alpha=c(NA,1,1,1,1,1,1,1,1,1,1,1,1),beta =c(0.5,0.1,1,1,1,1,1,1,1,1,1)) ```

Results