Понимание и исправление ошибок ложной сходимости в glmmTMB + lme4

Question

Я пытаюсь использовать glmmTMBчтобы запустить несколько итераций моей модели, но постоянно получать одну и ту же ошибку. Я попытался объяснить свой эксперимент ниже и вставил полную модель, которую пытаюсь запустить.

Предпосылки к эксперименту

Зависимая переменная, которую я пытаюсь смоделировать, - это количество копий бактериального гена 16S, которое в данном случае используется в качестве прокси для бактериальной биомассы.

Экспериментальный план состоит в том, что у меня есть отложения из 8 ручьев, которые падают по градиенту загрязнения (до нетронутого). (Фактор 1 = поток, с 8 уровнями).

Для каждого из 8 потоков было выполнено следующее: Осадок был добавлен в 6 лотков. 3 из этих лотков были помещены в канал искусственного потока, нагретого до 13°C, в то время как остальные 3 были нагреты до 17°C (фактор 2 = обработка с подогревом, с 2 уровнями). Всего имеется 16 каналов, и обработка Warming была назначена каналу случайным образом.

Затем я повторно измерил 3 лотка в каждом канале потока в четырех временных точках (Фактор 3 = День, с 4 уровнями).

В данный момент я рассматриваю лоток как настоящую биологическую копию, а не псевдореп, поскольку лотки расположены значительно дальше друг от друга в каналах, но это требует изучения.

Итак, чтобы подвести итог: условия модели (все указаны как факторы):

• Согревающая процедура (13 против 17oC)
• StreamID (1,2,3,4,5,6,7,8)
• День (T1, T4, T7, T14)

Полная модель, которую я предлагал, была такой:

X4_tmb.nb2<-glmmTMB(CopyNo~Treatment*Stream*Time, family=nbinom2, data=qPCR)   

Несмотря на то, что эта версия модели не включает случайный эффект, я хотел использовать glmmTMB пакет и не запускать его с помощью lme4, потому что я хотел изучить идею добавления компонента модели для учета дисперсии, а также изучить возможность добавления лотка как случайного эффекта (не уверен, что это правильно). Запустив все версии модели вglmmTMB, Я могу с уверенностью сравнить их результаты в AIC. Я бы не смог этого сделать, если бы запустил полную модель без компонента дисперсии вlme4 и другие, использующие glmmTMB.

К сожалению, для большинства итераций полной модели при использовании glmmTMB (я имею в виду последовательное удаление членов модели) я получаю такое же постоянное предупреждение:

Предупреждающее сообщение: In fitTMB(TMBStruc): проблема сходимости модели; ложная сходимость (8). См. Виньетку ("устранение неполадок")

Я попытался понять ошибку, но мне трудно понять, потому что сбивает с толку то, что когда я запускаю полную модель с помощью lme4, она работает без ошибок.

Это версия полной модели, которая работает в lme4,

X4 = glm.nb(CopyNo~Treatment*Stream*Time, data = qPCR

Насколько я понял из чтения https://www.biorxiv.org/content/10.1101/132753v1.full.pdf@ line 225, этот пакет можно использовать для перекрестного сравнения между GLM и GLMM. Вы знаете, правильно ли я это понял?

Я также использовал DHARMa пакет, чтобы помочь проверить модели и версию, которая не смогла сойтись с помощью glmmTMB, пройти тест KStest, тест дисперсии, тест выбросов и комбинированный тест скорректированного квантиля, но в идеале я не хочу ошибки сходимости.

Любая помощь будет принята с благодарностью.

Answer 1

Здесь куча.

предупреждение

К сожалению, с этим сложно что-то сделать: это общеизвестно неясное сообщение об ошибке. Как было предложено в Twitter, вы можете попробовать другой оптимизатор, например, включить

 control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS"))

в вашем звонке. Надеюсь, это даст очень похожий ответ (и в этом случае вы сделаете вывод, что предупреждение о конвергенции, вероятно, было ложным, так как разные оптимизаторы вряд ли потерпят неудачу одинаково) без предупреждения. (Вы можете попробоватьmethod="CG"выше в качестве третьей альтернативы.) (Обратите внимание, что есть небольшая ошибка с печатью и суммированием результатов при использовании альтернативных оптимизаторов, которая была недавно исправлена; вы можете установить версию для разработки, если вы работаете над этим до того, как исправление распространится на CRAN.)

модель "lme4"

В glm.nb()функция не изlme4 пакет, это из MASSпакет. Если бы у вас были случайные эффекты в модели, вы бы использовалиglmer.nb(), который находится вlme4 пакет... как в тестах переключения оптимизатора выше, если вы получите аналогичные ответы с glmmTMB а также glm.nb можно сделать вывод, что предупреждение от glmmTMB (собственно, это из nlminb() оптимизатор, который glmmTMB внутренние вызовы), вероятно, является ложным срабатыванием.

Самый простой способ проверить соизмеримость правдоподобия /AIC для разных пакетов - по возможности уместить одну и ту же модель в оба пакета, например

library(MASS)
library(glmmTMB)

quine.nb1 <- glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine)
quine.nb2 <- glmmTMB(Days ~ Sex/(Age + Eth*Lrn), data = quine,
                     family=nbinom2)
all.equal(AIC(quine.nb1),AIC(quine.nb2))  ## TRUE

другие детали

Одна из возможных проблем с вашей моделью заключается в том, что, подбирая полное трехстороннее взаимодействие трех категориальных переменных, вы пытаетесь оценить (2*4*8=) 64 параметра до 64*3=192 наблюдений (если я понимаю ваш экспериментальный дизайн правильно). Это и более вероятно приведет к численным проблемам (как указано выше) и может дать неточные результаты. Хотя некоторые рекомендуют это, я лично не рекомендую подход к выбору модели (все подмножества или последовательный, на основе AIC или на основе p-значения); Я бы предложил сделать Stream случайным эффектом, т.е.

CopyNo ~ Treatment + (Treatment|StreamID) + (1|Time/StreamID)

Это соответствует (1) общему лечебному эффекту; (2) изменение по потокам, изменение эффекта обработки по потокам; (3) изменение во времени и между потоками в определенные моменты времени. Здесь используется только 2 (фиксировано: перехват + обработка) + 3 (дисперсия перехвата, дисперсия обработки и их ковариация по потокам) + 2 (дисперсия во времени и между потоками во времени). Это не совсем "максимальная" модель; поскольку обе обработки измеряются каждый раз в каждом потоке, последним термином может быть (Обработка | Время / ID потока), но это значительно усложнит модель. Поскольку 4 группы - это немного для случайного эффекта, вы можете обнаружить, что хотите превратить последний член вTime + (1|Time:StreamID), который будет соответствовать времени как фиксированному и (потокам во времени) как случайному...