How to estimate multilevel Poisson models in lme4 with large count values?
<p><strong>The case:</strong>
I am trying to estimate gravity models of migration (<code>mig</code>) from 55 origin districts (<code>IDo</code>) to 54 destination districts (<code>IDd</code>). Two important predictors are the distance (<code>dist</code>) between each origin and destination district and an indicator variable of the contiguity (<code>contig</code>) of any origin-destination pair. The migration outcome variable is a measure of the count of people migrating between pairs of districts (ranging from zero to a couple of thousands). Below is a test data set that approximates my data situation. In the real data, the migration outcome variable is slightly less ideally poisson distributed (<code>min=0, max=9450, mean=85, median=10</code>).</p>

<pre><code>library(lme4)

#*** Generate test data set
set.seed=777
td=data.frame(IDo=rep(1:55,each=55),IDd=rep(1:55,times=55),dist=runif(3025,0.186,12.7),contig=rbinom(3025,1,p=0.08), stringsAsFactors=F)
td=td[td$IDo!=td$IDd,]              # remove cases for which origin and destination are the same
td=rbind(td,td)                     # dublicate data for two years
td$year=rep(c(1,0),each=2970)       # generate year dummy variable
td$mig=rpois(5940, lambda = 1)*1000 # generate migrant count variable

# Cross-classified random effects Poisson models
m1=formula(mig~dist+contig+year+(1|IDo)+(1|IDd))
fm0=glmer(m1,data=td,family="poisson",control=glmerControl(optimizer="bobyqa"), nAGQ=0) # Adaptive Gauss-Hermite Quadrature
fm1=glmer(m1,data=td,family="poisson",control=glmerControl(optimizer="bobyqa"), nAGQ=1) # Laplace Approximation

# Regular Poisson models
m2=formula(mig~dist+contig+year+factor(IDo)+factor(IDd))
fm2=glm(m2,data=td,family="poisson")
</code></pre>

<p><strong>The problem:</strong>
I am running cross-classified multilevel Poisson models using <code>lme4</code>. When using the default nAGQ=1 setting (<code>fm1</code>), I always get the following warning message, and in the real data case, the model fails to converge (<code>In checkConv(attr(opt, "derivs"), opt$par, control$checkConv Model failed to converge with max|grad| = 0.00248588 (tol = 0.001, component 1)</code>).</p>

<pre><code>Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?
</code></pre>

<p>However, this problem does not emerge with setting of nAGQ=0 (<code>fm0</code>) or when running it as a regular fixed effects model (<code>fm2</code>).</p>

<p><em>Does anyone have a suggestion what to do to avoid the problem that underlies this warning message (and leads to convergence issues in the real data)?</em> The predictor variables have all a reasonably small scale, so no need to rescale these variables. I could potentially rescale my outcome variable to have smaller values (<code>td$mig=round(td$mig/1000,0)</code>) but this would change the distribution and should not be done as indicated in <a href="https://stats.stackexchange.com/questions/270/poisson-regression-with-large-data-is-it-wrong-to-change-the-unit-of-measuremen">this post</a>.</p>