如何估计 lme4 中具有大计数值的多级泊松模型？

Question

案例： 我正在尝试估计从 55 个起源区 (IDo) 到 54 个目的地区 (IDd) 的迁移引力模型 (mig)。两个重要的预测变量是每个起点和终点区域之间的距离 (dist) 和任何起点-终点对的连续性 (contig) 的指示变量。迁移结果变量是衡量在成对地区之间迁移的人数（范围从零到几千）。下面是一个近似于我的数据情况的测试数据集。在真实数据中，迁移结果变量的泊松分布稍差（min=0, max=9450, mean=85, median=10）。

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")

问题： 我是运行使用 lme4 的交叉分类多级泊松模型。当使用默认的 nAGQ=1 设置（fm1）时，我总是收到以下警告消息，而在真实数据情况下，模型无法收敛（In checkConv(attr(opt, "derivs"), opt$par, control$checkConv Model failed to converge with max|grad| = 0.00248588 (tol = 0.001, component 1)）。

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

但是，设置 nAGQ=0 (fm0) 或运行作为常规固定效应模型 (fm2) 时不会出现此问题。

有没有人建议如何避免此警告消息背后的问题（并导致实际数据中的收敛问题）？预测变量具有所有规模相当小，因此无需重新调整这些变量。我可能会重新调整我的结果变量以具有更小的值 (td$mig=round(td$mig/1000,0))，但这会改变分布，不应按照 this post.

中的指示进行

Answer 1

我真的认为这是误报；拥有如此大量的计数数据是不寻常的……您是否尝试过使用不同的优化器进行拟合（参见 ?lme4::convergence）？（我刚刚用 control=glmerControl(optim="nloptwrap") 尝试了你的例子，得到了几乎相同的结果。）

我有点担心你的分布假设。我知道这只是一个可重现的例子，但是绘制一个泊松样本并将其乘以 1000 确实 not 给你一个泊松分布的值......这是一个模拟分布的小后验预测模拟拟合模型下响应的第 90 个分位数 - way 不同于观察值 ...

ss <- simulate(fm1,1000,seed=101)
qq <- sapply(ss,quantile,0.9)
hist(qq,breaks=50,col="gray")
summary(qq)
##   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1106    1153    1167    1168    1182    1255 
quantile(td$mig,0.9) ## 2000

如何估计 lme4 中具有大计数值的多级泊松模型？

How to estimate multilevel Poisson models in lme4 with large count values?

r

lme4

poisson