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:
The age group parameters blow out of proportion starting with age group 8. Why is this happening? Where did I go wrong?
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