Как подобрать смешанную модель с переменной отклика от 0 до 1?

Я пытаюсь использовать lme4::glmer()для подгонки биномиальной обобщенной смешанной модели (GLMM) с зависимой переменной, которая является не двоичной, а непрерывной переменной от нуля до единицы. Можно думать об этой переменной как о вероятности; на самом деле это вероятность того, как сообщили человеческих субъектов (в эксперименте , который я помочь анализирующего). Т.е. это не «дискретная» дробь, а непрерывная переменная.

Мой glmer()звонок не работает должным образом (см. Ниже). Почему? Что я могу сделать?

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

Подробнее

По-видимому, можно использовать логистическую регрессию не только для двоичного DV, но также и для непрерывного DV от нуля до единицы. Действительно, когда я бегу

glm(reportedProbability ~ a + b + c, myData, family="binomial")

Я получаю предупреждение

Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!

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

Однако то, что я на самом деле хочу использовать,

glmer(reportedProbability ~ a + b + c + (1 | subject), myData, family="binomial")

Это дает мне то же самое предупреждение, возвращает модель, но эта модель явно не в порядке; оценки фиксированных эффектов очень далеки от оценок glm()и от средних по предмету. (И мне нужно включить glmerControl(optimizer="bobyqa")в glmerвызов, иначе он вообще не сходится.)

Имеет смысл начать с более простого случая отсутствия случайных эффектов.

Существует четыре способа справиться с непрерывной переменной ответа «ноль к одному», которая ведет себя как дробь или вероятность ( это наша самая каноническая / просматриваемая / просматриваемая ветка по этой теме, но, к сожалению, здесь не все четыре варианта обсуждаются):

1. $p=m/n$$n$n$N$

   glm(p ~ a+b+c, myData, family="binomial", weights=n)

2. $p$$p$$0$$1$

   betareg(p ~ a+b+c, myData)

3. Логите преобразование ответа и используйте линейную регрессию. Обычно это не рекомендуется.

   lm(log(p/(1-p)) ~ a+b+c, myData)

4. Подберите биномиальную модель, но затем рассчитайте стандартные ошибки с учетом чрезмерной дисперсии. Стандартные ошибки могут быть вычислены различными способами:

   • (а) масштабированные стандартные ошибки с помощью оценки избыточной дисперсии ( один , два ). Это называется "квазибиномиальным" GLM.

   • (б) устойчивые стандартные ошибки с помощью сэндвич-оценки ( один , два , три , четыре ). Это называется «дробный логит» в эконометрике.

   
   (A) и (b) не являются идентичными (см. Этот комментарий и разделы 3.4.1 и 3.4.2 в этой книге , и этот пост SO, а также этот и этот ), но, как правило, дают схожие результаты. Вариант (а) реализован glmследующим образом:

   glm(p ~ a+b+c, myData, family="quasibinomial")

Те же четыре способа доступны со случайными эффектами.

1. Используя weightsаргумент ( один , два ):

   glmer(p ~ a+b+c + (1|subject), myData, family="binomial", weights=n)

   Согласно второй ссылке выше, это может быть хорошей идеей для моделирования избыточной дисперсии, см. Там (а также # 4 ниже).

2. Использование бета-смешанной модели:

   glmmadmb(p ~ a+b+c + (1|subject), myData, family="beta")

   или

   glmmTMB(p ~ a+b+c + (1|subject), myData, 
           family=list(family="beta",link="logit"))

   Если в ответных данных есть точные нули или единицы, то можно использовать бета-модель с нулевой / одной раздувкой в glmmTMB.

3. Используя логит-преобразование ответа:

   lmer(log(p/(1-p)) ~ a+b+c + (1|subject), myData)

4. Учет избыточной дисперсии в биномиальной модели. Это использует другой прием: добавление случайного эффекта для каждой точки данных:

   myData$rowid = as.factor(1:nrow(myData))
   glmer(p ~ a+b+c + (1|subject) + (1|rowid), myData, family="binomial",
         glmerControl(optimizer="bobyqa"))

   По некоторым причинам это не работает должным образом, так как glmer()жалуется на нецелые числа pи дает бессмысленные оценки. Решение, которое я придумал, состоит в том, чтобы использовать поддельную константу weights=kи убедиться, что p*kона всегда целочисленная. Это требует округления, pно выбор kдостаточно большого значения не должен иметь большого значения. Результаты, похоже, не зависят от стоимости k.

   k = 100
   glmer(round(p*k)/k ~ a+b+c + (1|subject) + (1|rowid), myData, 
         family="binomial", weights=rowid*0+k, glmerControl(optimizer="bobyqa"))

   Позднее обновление (январь 2018 г.): это может быть неверный подход. Смотрите обсуждение здесь . Я должен исследовать это больше.

В моем конкретном случае вариант № 1 недоступен.

Вариант № 2 очень медленный и имеет проблемы с конвергенцией: glmmadmbдля его запуска требуется пять-десять минут (и он все еще жалуется, что он не сходится!), Тогда как он lmerработает за доли секунды и glmerзанимает пару секунд. Обновление: я попробовал, glmmTMBкак предложено в комментариях @BenBolker, и он работает почти так же быстро, как glmerи без проблем сходимости. Так что это то, что я буду использовать.

Варианты № 3 и № 4 дают очень похожие оценки и очень похожие доверительные интервалы Вальда (полученные с confint ). Я не большой поклонник # 3, хотя, потому что это своего рода обман. И № 4 чувствует себя немного хакером.

Огромное спасибо @Aaron, который указал мне на №3 и №4 в своем комментарии.