Понимание алгоритма блеска;помогите с отладкой используя R

Question

Я пытаюсь скопировать некоторое подмножество lme4 в Haskell. Я справился с частью lmer, более или менее, с результатами, которые соответствуют примерам Бейтса (2009) «Реализация линейной смешанной модели в lme4». Так что теперь я перехожу к обобщенному. Немного обнадеживает то, что я могу выполнить свой блестящий эквивалент на примерах lmer, используя семейство Гаусса и функцию тождественной ссылки. Так что все не 10000 * дико неправильно. Но когда я пытаюсь соответствовать примеру cbpp, даже очень упрощенная версия, где я просто пытаюсь соответствовать

example <- glmer(incidence/size ~ 1 + (1 | herd), data = cbpp,
                 family = binomial(), weights = size)

Все не работает. Под этим я подразумеваю, что моя версия на Haskell (https://github.com/adamConnerSax/glm-haskell/blob/GLMMs/test/glm/GLMM-test.hs) не сходится к одному и тому же ответу.

Было бы очень полезно две вещи:

1. Isесть лучшая ссылка на алгоритм glmer, чем в различных статьях Бейтса (2007, 2009, 2018)? Все они имеют много специфических особенностей в части LMM, но намного меньше в бите GLMM.

2. Предположим, я выполнил приведенный выше код с результатом в example. В lme4, есть ли способ просто получить значение различных функций отклонения для данных beta и u вместо решения? Могу ли япросто запустите итеративно взвешенные наименьшие квадраты для фиксированного значения theta, вместо того, чтобы запускать весь цикл оптимизации? Чем больше у меня точек сравнения между R и моим кодом на Haskell, тем проще для меня это будетчтобы увидеть, где я ошибся.

Answer 1

1. Вы можете получить незаконченный черновик с подробной теорией и реализацией GLMM с пакетом lme4 на r-forge (большинство lme4 разработок перенесено с r-подделать GitHub, но черновик этого документа все еще здесь)

2. Вы можете абсолютно просто получить функцию отклонения. Самый простой способ - это взять подобранную модель и запустить update(fitted_model, devFunOnly=TRUE) (или вы можете взять тот же набор аргументов и добавить devFunOnly=TRUE в первую очередь). Есть несколько потенциальных хитрых битов о том, выполняется ли первый (nAGQ=0) или второй проход (с nAGQ>0);см. ?modular.

   Однако эта функция отклонения может быть не совсем такой, как вы хотите, поскольку она принимает beta и theta - но не u - в качестве ввода. В рамках функции отклонения функция pwrssUpdate запускает алгоритм PIRLS. Но это, по крайней мере, приблизит вас ... вы можете посмотреть на исходный код или, если вы сохранили функцию отклонения devfun, вы можете посмотреть на environment(devfun)$pwrssUpdate