данные о сверхдисперсии с нулевым накачиванием, ошибка glmmTMB в R

Question

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

Перед тем как работать с данными, я проверил их на нормальность (они не завышены), однородность дисперсии, корреляции и выбросы. У данных было два выброса, которые я удалил из набора данных linekd выше. Есть 351 наблюдение из 18 точек (prop_id).

Данные выглядят так:

euc0 ea_grass ep_grass np_grass np_other_grass month year precip season   prop_id quad
 3      5.7      0.0     16.7            4.0     7 2006    526 Winter    Barlow    1
 0      6.7      0.0     28.3            0.0     7 2006    525 Winter    Barlow    2
 0      2.3      0.0      3.3            0.0     7 2006    524 Winter    Barlow    3
 0      1.7      0.0     13.3            0.0     7 2006    845 Winter    Blaber    4
 0      5.7      0.0     45.0            0.0     7 2006    817 Winter    Blaber    5
 0     11.7      1.7     46.7            0.0     7 2006    607 Winter    DClark    3

Переменная ответа euc0 и случайные эффекты prop_id а также quad_id. Остальные переменные являются фиксированными эффектами (все они представляют собой процентное покрытие различных видов растений).

Модель, которую я хочу запустить:

library(glmmTMB)
seed0<-glmmTMB(euc0 ~ ea_grass + ep_grass + np_grass + np_other_grass + month + year*precip + season*precip + (1|prop_id)  + (1|quad), data = euc, family=poisson(link=identity))

fit_zinbinom <- update(seed0,family=nbinom2) #allow variance increases quadratically

Ошибка, которую я получаю после запуска seed0 код:

Ошибка в optimHess(par.fixed, obj$fn, obj$gr): градиент в optim оценивается как длина 1, а не 15. Вдобавок: было 50 или более предупреждений (используйте warnings(), чтобы увидеть первые 50)

warnings() дает:

1. In (function (start, objective, gradient = NULL, hessian = NULL,  ... :
  NA/NaN function evaluation

Я также обычно имею в виду центр и стандартизацию моих числовых переменных, но это устраняет только первую ошибку и сохраняет NA/NaNошибка. Я пробовал добавитьglmmTMBControlтакой оператор, как этот OP, но он просто открыл целый новый мир ошибок.

Как я могу это исправить? Что я делаю не так?

Будем признательны за подробное объяснение, чтобы я мог научиться лучше решать эту проблему в будущем. В качестве альтернативы я открыт дляMCMCglmm решение, поскольку эта функция также может работать с данными такого рода (несмотря на то, что на выполнение требуется больше времени).

Answer 1

Неполный ответ...

• модели связи идентичности для распределений ответов в ограниченной области (например, гамма или Пуассон, где отрицательные значения невозможны) являются проблематичными с вычислительной точки зрения; на мой взгляд, они часто концептуально проблематичны, хотя в их пользу есть несколько разумных аргументов. У вас есть для этого веская причина?
• Это довольно небольшой набор данных для модели, которую вы пытаетесь подогнать: 13 предикторов с фиксированным эффектом и 2 предиктора со случайным эффектом. Эмпирическое правило состоит в том, что вам нужно примерно в 10-20 раз больше наблюдений: это, кажется, согласуется с вашими 345 или около того наблюдениями, но... только 40 ваших наблюдений ненулевые! Это означает, что ваше "эффективное" количество наблюдений / количество информации будет намного меньше (см. Стратегии регрессионного моделирования Фрэнка Харрелла для более подробного обсуждения этого вопроса).

Тем не менее, позвольте мне рассказать о некоторых вещах, которые я пробовал, и о том, чем я закончил.

• GGally::ggpairs(euc, columns=2:10) не обнаруживает ничего явно ужасного в данных (я выбросил точку данных с euc0==78)

Чтобы попытаться заставить работать модель идентификации-связи, я добавил код в glmmTMB. У вас должна быть возможность установить черезremotes::install_github("glmmTMB/glmmTMB/glmmTMB@clamp")(обратите внимание, что для установки вам потребуются компиляторы и т. д.). Эта версия принимает отрицательные предсказанные значения и заставляет их быть неотрицательными, добавляя при этом соответствующий штраф к отрицательной логарифмической вероятности.

Используя новую версию glmmTMB, я не получаю сообщения об ошибке, но получаю следующие предупреждения:

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

Неположительно определенная матрица Гессе (вторая производная) означает, что существуют некоторые (все еще трудно устраняемые) проблемы. heatmap(vcov(f2)$cond,Rowv=NA,Colv=NA)Позвольте мне взглянуть на ковариационную матрицу. (Мне также нравитсяcorrplot::corrplot.mixed(cov2cor(vcov(f2)$cond),"ellipse","number"), но это не работает, когда vcov(.)$condнеположительно определен. В крайнем случае вы можете использоватьsfsmisc::posdefify() чтобы заставить его быть положительно определенным...)

Пробовал масштабировать:

eucsc <- dplyr::mutate_at(euc1,dplyr::vars(c(ea_grass:precip)), ~c(scale(.)))

Это поможет некоторым - прямо сейчас мы все еще делаем несколько глупых вещей, например, рассматриваем год как числовую переменную, не центрируя ее (так что "перехват" модели находится в году 0 по григорианскому календарю...)

Но это все еще не решает проблему.

Если присмотреться к ggpairs сюжет, похоже season а также year сбиты с толку: with(eucsc,table(season,year)) показывает, что наблюдения происходят весной и зимой в один год и осенью в другой год. season а также month также сбиты с толку: если мы знаем месяц, мы автоматически знаем время года.

На этом этапе я решил отказаться от идентификационной ссылки и посмотреть, что произошло. update(<previous_model>, family=poisson)(т.е. использование Пуассона со стандартной ссылкой журнала) сработало! Так было с использованиемfamily=nbinom2, что было намного лучше.

Я просмотрел результаты и обнаружил, что доверительные интервалы для коэффициентов сезона осадков X были сумасшедшими, поэтому отказался от члена взаимодействия (update(f2S_noyr_logNB, . ~ . - precip:season)), после чего результаты выглядят разумными.

Несколько заключительных замечаний:

• дисперсия, связанная с квадратом, фактически равна нулю
• Я не думаю, что вам обязательно нужна нулевая инфляция; низкие средние и чрезмерная дисперсия (т.е.family=nbinom2), вероятно, достаточно.
• распределение остатков выглядит нормально, но, похоже, некоторые модели не подходят (library(DHARMa); plot(simulateResiduals(f2S_noyr_logNB2))). Я бы потратил некоторое время на построение графиков остатков и прогнозируемых значений против различных комбинаций предикторов, чтобы увидеть, сможете ли вы локализовать проблему.

PS Более быстрый способ увидеть, что с фиксированными эффектами что-то не так (мультиколлинеарность):

X <- model.matrix(~ ea_grass + ep_grass +
                   np_grass + np_other_grass + month +
                   year*precip + season*precip,
                  data=euc)
ncol(X)  ## 13
Matrix::rankMatrix(X) ## 11

lme4 есть такие тесты и механизмы для автоматического удаления столбцов с псевдонимами, но они не реализованы в glmmTMB в настоящий момент.