данные о сверхдисперсии с нулевым накачиванием, ошибка glmmTMB в R
Я работаю с данными подсчета (доступны здесь), которые являются нулевыми и чрезмерно диспергированными и имеют случайные эффекты. Пакет, который лучше всего подходит для работы с такими данными, - это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
решение, поскольку эта функция также может работать с данными такого рода (несмотря на то, что на выполнение требуется больше времени).
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
в настоящий момент.