Как оценить многоуровневые модели Пуассона в lme4 с большими значениями счета?

Question

Случай: пытаюсь оценить гравитационные модели миграции (mig) из 55 районов происхождения (IDo) до 54 районов назначения (IDd). Двумя важными предикторами являются расстояние (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)) но это изменит распределение и не должно быть сделано, как указано в этом посте.

Answer 1

Я действительно думаю, что это ложный позитив; необычно иметь данные подсчета с таким большим количеством... пробовали ли вы подгонять к различным оптимизаторам (см. ?lme4::convergence) (Я только что попробовал ваш пример с control=glmerControl(optim="nloptwrap") и получил почти идентичные результаты.)

Я немного беспокоюсь о ваших предположениях о распределении. Я знаю, что это всего лишь воспроизводимый пример, но рисование образца Пуассона и умножение его на 1000 не дает вам значение, распределенное по Пуассону... вот небольшая апостериорная прогностическая симуляция, которая имитирует распределение 90-го квантиля ответа под модель - которая сильно отличается от наблюдаемого значения...

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