Using mle2() for age-period-cohort models

Question

I've followed guidelines for comparing models in Chapter 6 of Bolker's Ecological Models and Data in R, applying code used in this section to cancer count data. The models include parameters for age group (a) and birth cohort (c), and are as follows:

# Value of lambda is the same for all age groups and birth cohorts
poisfit.0 <- mle2(counts ~ dpois(lambda = a*c*pyrs), start=list(a=3,c=0.7),data=x)  

Warning message: In dpois(x, lambda, log) : NaNs produced

Coefficients: Estimate Std. Error z value Pr(z)
a 3.1932757 0.0022079 1446.28 < 2.2e-16
c 0.8429316 0.0083643 100.78 < 2.2e-16

-2 log L: 17550.32

# Use coefficients from simplest model to build more complex models  
start.ab <- coef(poisfit.0)

# Model with varying parameter values for age groups and birth cohorts
poisfit.ab = mle2(counts ~ dpois(lambda = a*c*pyrs), start = start.ab, data = x, parameters=list(a~x\$ageGroup, c~x\$bCohort))

# Same as above, but without intercept
poisfit.ab1 = mle2(counts ~ dpois(lambda = a*c*pyrs), start = start.ab, data = x, parameters=list(a~x\$ageGroup-1, c~x$bCohort-1))

The data look like this, and go on for 110 rows (14 age groups, 16 birth cohorts):
ageGroup bCohort pyrs counts
1 13 1.4006 0
1 14 32.6925 1
1 15 49.7632 1
1 16 49.7059 0
2 12 1.2528 0
2 13 30.3879 3
2 14 47.7709 1
2 15 52.3699 0
2 16 55.7669 2
3 11 1.0257 0

A sample of coefficients for the model lacking intercept:

Estimate    Std.    Error   z   value   Pr(z)   
a.x\$ageGroup1  0.101115    0.072032    1.4037  0.1603938   
a.x\$ageGroup2  0.189811    0.078609    2.4146  0.0157517       
a.x\$ageGroup3  0.335937    0.086470    3.8850  0.0001023   
a.x\$ageGroup4  0.760921    0.117493    6.4763  9.400e-11   
a.x\$ageGroup5  1.282271    0.141785    9.0438  <   2.2e-16
a.x\$ageGroup6  2.221877    0.180402    12.3162 <   2.2e-16
a.x\$ageGroup7  3.694786    0.227830    16.2173 <   2.2e-16
a.x\$ageGroup8  6.629915    0.313676    21.1362 <   2.2e-16
a.x\$ageGroup9  10.103253   0.400364    25.2351 <   2.2e-16
a.x\$ageGroup10 16.198128   0.537537    30.1340 <   2.2e-16
a.x\$ageGroup11 24.141987   0.686824    35.1502 <   2.2e-16
a.x\$ageGroup12 32.166317   0.811856    39.6207 <   2.2e-16
a.x\$ageGroup13 39.877367   0.941045    42.3756 <   2.2e-16
a.x\$ageGroup14 36.982237   1.169403    31.6249 <   2.2e-16
c.x\$bCohort1   0.058183    0.026073    2.2315  0.0256476   
c.x\$bCohort2   0.135923    0.023366    5.8171  5.988e-09   

There are two questions I have about this model:

1. The age group parameters blow out of proportion starting with age group 8. Why is this happening? Where did I go wrong?

2. The starting parameter values in the simple model poisfit.0 are both positive. They must both have the same sign in order for mle2() to accept them. If one is positive and the other negative, mle2() issues this error:

Error in optim(par = c(-0.2, 3.4), fn = function (p) : initial value in 'vmmin' is not finite

I don't know enough about what goes on under the hood to understand why this happens. I'd appreciate insight into these questions.

Thank you,

Susana

Answer 1

There are several important problems here.

• Are you sure you can't just use Poisson regression (i.e. glm(counts~ageGroup+bCohort,data=x,family=poisson) (i.e. see ch 9 of Bolker 2008, or any other book on GLMs, such as Faraway 2006)?
• In your first model, a and c are jointly unidentifiable -- for example, a=1, c=2 would have exactly the same likelihood as a=2, c=1 (or any other combination with a*c=2)
• you should probably use (a~ageGroup c~bCohort) instead of including the x$
• I don't think the parameters are "blowing up", it's just that there is a progressive trend through age (you should seriously consider a model that treats age group as a continuous predictor)
• the problem you refer to with the starting values is that you can't start with parameters that make the lambda value negative. (mle2, or rather most of the underlying optimizers in optim, can handle parameters that return NA for the negative log-likelihood during the course of the optimization, but not right at the beginning of the optimization). You could fit lambda on the log scale (i.e. loglambda=loga+logc), or use bounded optimization.