Нетехническое продолжение LMM

Итак, недавно мы делали нетехническое введение в LMM, а значит, пришло время продолжить этот разговор.

Нужные нам пакеты: pandas, statsmodels

Ссылка на данные: тык

Приступим!

Больше случайных эффектов!

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

Перекрывающиеся эффекты касаются тех случаев, когда они [эффекты] независимы друг от друга. Если представить, что есть один случайный эффект с подгруппами A, B, C и другой с подгруппами 1, 2, 3, то их взаимодействие можно изобразить так:

Как видите, эти подгруппы могут встречаться в любых комбинациях, и они всегда одинаковы. Конечно, вариант выше является идеальным и чаще всего встречается в сбалансированных экспериментальных исследованиях, но в реальной жизни случайные эффекты перекрываются лишь частично (например, C встречается только в комбинациях с 1 и 3, а B только с 2).

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

Визуально это выглядит подобным образом:

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

Не дайте себя обмануть! Подгруппа 1, вложенная в подгруппу А, и подгруппа 1, вложенная в подгруппу B – это две разные подгруппы и их нельзя путать. Обычно статистические пакеты могут сами их различать при правильной спецификации модели, что всё ещё требует понимания данных с вашей стороны, но об этом дальше.

Что по нашим данным?

В нашем датасете именно такая ситуация: помимо разделения на школы, ученики разделяются ещё и по классам. Это классический пример вложенных эффектов - класс 1А из школы №1 и класс 1А из школы №2 отличаются друг от друга и состоят из разных людей!

С такой ситуацией statsmodels справляется без проблем благодаря аргументу vc_formula – на вход он принимает словарь вида {‘название_эффекта_в_модели_1’ : ‘0 + название_эффекта_в_данных_1’, ‘название_эффекта_в_модели_2’ : ‘0 + название_эффекта_в_данных_2’}. Обратите внимание на 0 – так модель поймёт, что ей не нужно рассчитывать ещё один свободный член. Также нужно указать аргумент re_formula = '1', чтобы строился случайный свободный член, связанный со всей нашей структурой:

vc = {'classid': '0 + C(classid)'} #C() для указания категориального характера данных
model_nested = smf.mixedlm('mathgain ~ mathkind + ses + mathprep',
             data = classes, groups = classes['schoolid'], 
             vc_formula = vc, re_formula = '1').fit(method = 'lbfgs')
             
model_nested.summary()

В качестве результата получится следующее:

Появился дополнительный параметр classid Var, который отражает дисперсию свободных членов по классам. Как можно заметить, она почти такая же, как и по школам! Хотя нельзя не заметить, что дисперсия свободных членов по школам сильно ниже, чем в прошлой модели - видимо, изначально школы «присвоили» себе часть вариабельности, которая относится к классам.

Построим её не-REML версию и сравним с невложенной моделью:

model_nested_2 = smf.mixedlm('mathgain ~ mathkind + ses + mathprep',
               data = classes, groups = classes['schoolid'],
               vc_formula = vc, re_formula = '1').fit(method = 'lbfgs', reml = False)

И мы можем убедиться, что вложенная модель действительно лучше модели с одним случайным эффектом!

Уникальность кодов

Напоминаем, что эффекты низшего порядка могут быть как уникальными, так и неуникальными – например, 1 класс может повториться сразу в нескольких школах. Если эффектов только два, то statsmodels разбирается с этой проблемой автоматически. Мы можем дополнительно убедиться в этом, создав дополнительную переменную из комбинаций id школ и id классов:

classes['class_unique'] = classes.classid.astype('str') + '/' +  classes.schoolid.astype('str')

В результате у нас получается столбец с уникальными значениями для каждого сочетания класс/школа:

Построим модель с новой переменной:

vc = {'class_unique': '0 + C(class_unique)'}
model_nested_3 = smf.mixedlm('mathgain ~ mathkind + ses + mathprep',
               data = classes, groups = classes['schoolid'],
               vc_formula = vc, re_formula = '1').fit(method = 'lbfgs')
               
model_nested_3.summary()

Результат идентичен предыдущему:

Однако если вы захотите использовать большее количество эффектов, каждый из которых вложен в другой, то уникальность кодов внутри каждого уровня иерархии обязательна!

А что с перекрывающимися?

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

Наши данные были бы перекрывающимися только в том случае, если 1А класс (и родственные ему) был единым многомерным созданием, одновременно существующим в разных школах и состоящим из одинаковых детей. Тогда, правда, более интересными были бы вопросы «разумно ли оно?» и «что оно хочет от наших детей?», но уже за пределами данной статьи. А код для перекрывающейся модели был бы такой:

classes['group'] = 1 #создаем переменную из одних единиц

vc = {'classid': '0 + C(classid)', 'schoolid': '0 + C(schoolid)'}
model_crossed = smf.mixedlm('mathgain ~ mathkind + ses + mathprep',
              data = classes, groups = classes['group'],
              vc_formula = vc, re_formula = '1').fit(method = 'lbfgs')

GEE: a challenger appears

Возможно, вы сейчас всё это читаете и думаете «а зачем же мне все эти сложности со случайными эффектами, когда я хочу знать только фиксированные?». Если это действительно так, то вам может понравиться альтернативный метод моделирования сгруппированных данных под названием Generalized Estimating Equations (GEE). В чём же разница?

Метод не концентрируется на индивидуальной вариабельности и даёт только «фиксированные» эффекты. Таким образом, мы имеем более простую и однозначную модель.
Метод полупараметрический: хотя он всё ещё опирается на допущения используемых распределений, ряд внутренних расчётов не связан ни с каким распределением. В результате получается модель, более устойчивая к нарушениям своих допущений по сравнению с LMM.
Метод подходит только для одной разбивки по подгруппам...в большинстве случаев. Позднее мы расскажем про исключение.
Метод позволяет выбирать внутригрупповую ковариационную структуру. На этом остановимся поподробнее.

Ковариационные структуры

Вариационно-ковариационные матрицы являются важной составной частью любой линейной модели – они задают характер, нестандартизованную силу взаимодействия между переменными модели, их дисперсию, а в контексте моделирования сгруппированных данных касаются важного вопроса: насколько похожи друг на друга представители одной подгруппы?

Самый простой, кажущийся логичным, вариант – все коррелируют/ковариируют по-разному. В конце концов, люди разные, вариабельности у них много, кто-то похож друг на друга меньше, кто-то больше… Эта идея лежит в основе неограниченной структуры – у всех разные ковариации. При всей теоретической привлекательности такой идеи мы натыкаемся на ужасную практику: количество параметров в таких моделях растёт очень быстро. В результате может получиться так, что наша модель построится примерно никогда – параметров попросту настолько много, что никаких мощностей компьютера не хватит на расчёты.

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

Промежуточный вариант – взаимозаменяемая структура, которую можно увидеть под названием сложной симметрии. Тут представители одной группы коррелируют/ковариируют между собой одинаково – один коэффициент на всех.

Есть и множество других ковариационных структур - например, ауторегрессивная для моделирования изменяющихся во времени данных, разные варианты пространственных для моделирования географических данных и так далее. GEE позволяют свободно выбирать между разными вариантами этих структур, в отличие от LMM, которые обычно завязаны только на одной. Хотя в некоторых других статпакетах, не связанных с Python, можно менять их и в LMM.

Как это сделать в Python?

Это выполняется через общую функцию gee(), а разные виды ковариационных структур можно задавать через подмодуль sm.cov_struct с подачей результата в аргумент cov_struct. В остальном использование практически идентично коду модели случайного свободного члена, как мы и увидим далее.

Независимая структура

Вернёмся к нашим школам. Для начала попробуем смоделировать наши данные как независимые внутри школ.

ind = sm.cov_struct.Independence() #задаём независимую структуру
gee_ind = smf.gee('mathgain ~ mathkind + ses + mathprep',
        data = classes, groups = classes['schoolid'],
        cov_struct = ind).fit()

gee_ind.summary()

Посмотрим на результат:

Можем сравнить результат с классической линейной моделью, которую мы подсчитали в самом начале. Как видите, коэффициенты идентичны, но остальные параметры иные:

Вдобавок к этому, мы можем использовать метод .summary() на объекте с самой ковариационной структурой! Результат довольно предсказуем:

Взаимозаменяемая структура

Попробуем другую структуру. Теперь смоделируем эти наблюдения как одинаково скоррелированные!

ex = sm.cov_struct.Exchangeable()
gee_ex = smf.gee('mathgain ~ mathkind + ses + mathprep',
       data = classes, groups = classes['schoolid'],
       cov_struct = ex).fit()
       
gee_ex.summary()

Результаты изменились, но не то чтобы сильно:

Зато метод объекта с ковариационной структурой говорит более интересные вещи:

Согласно этой модели, ученики в рамках одной школы не так уж и похожи друг на друга! Видимо, поэтому коэффициенты и не изменились так сильно.

Вложенная структура

Помните, мы говорили про то, что GEE работают только с одной разбивкой по подгруппам? Настало время того самого исключения! Эта ковариационная структура позволяет создать подобие LMM со вложенными эффектами (перекрывающиеся всё ещё недоступны). Также появляется дополнительный аргумент dep_data - он принимает на вход строку со всеми подгруппами через «+» в порядке «от высшего к низшему». Здесь тоже важно поставить в начале 0, чтобы модель не строила лишний свободный член:

nest = sm.cov_struct.Nested()
gee_nest = smf.gee('mathgain ~ mathkind + ses + mathprep',
         data = classes, groups = classes['schoolid'],
         cov_struct = nest, dep_data = '0 + classid').fit()
         
gee_nest.summary()

Но получаем мы всё ещё только фиксированные эффекты:

Зато смотрите, о чём нам говорит структура:

Мы видим внутригрупповую дисперсию на каждом уровне. Дисперсия по классам немного ниже, чем дисперсия по школам – значит, классы всё же более гомогенны, и ученики в них больше похожи друг на друга! Обратите внимание на огромную остаточную (residual) дисперсию – огромное количество вариабельности в данных так и не было объяснено, и нужно искать новые важные предикторы.

Какая модель лучше?

Вы не поверите, но мы снова не можем пользоваться AIC! В данном случае ситуация ещё хуже: из-за полупараметрического характера модели она использует т.н. функцию псевдоправдоподобия, которая для классических информационных критериев не подходит в принципе.

Благо для таких моделей был создан псевдоинформационный критерий (quasi-information criterion, QIC) со схожей интерпретацией. Более того, statsmodels даёт его в двух версиях: классической и упрощённой (QICu), которая не может использоваться для сравнения ковариационных структур, но подходит для различий в других параметрах модели. А вы думали, это у вас синдром самозванца :)

Реализуется он с помощью метода .qic() с дополнительным фокусом: нужно также указать параметр scale для той модели, у которой больше всего коэффициентов. Если этого не делать, то он использует для каждой модели свой параметр и пообещает страшные кары. Чтобы этого не произошло, возьмём scale у последней модели:

Нормальный QIC здесь выражен первым числом, но особой разницы нет: самой лучшей моделью оказалась независимая. Обратите внимание, что этот вывод расходится с выводом по LMM, где лучшей оказалась как раз вложенная. Это может быть связано как с разными допущениями и алгоритмами двух классов моделей, так и с особенностями самого QIC. «Правильный» ответ зависит от того, что вас больше интересует в данных, но это именно то решение, которое мы приняли при выборе изначального метода анализа. Других вариантов понять «истину» нет. Совсем нет. Честно.

При чём тут GLM?

Одна из приятных особенностей обоих названных классов моделей заключается в том, что существуют их GLM-версии! Для тех, кто не знает: обобщённые линейные модели (generalized linear models, GLM) работают практически как обычные, но вместо нормального распределения предполагают какое-то другое (биномиальное, Пуассона, Вейбулла и т.д.). Не так давно был про это вебинар, так что рекомендуем его посмотреть.

Обрадовались? А теперь время горькой, печальной правды: в statsmodels GLM и LMM не очень хорошо дружат. Эти модели там существуют, но есть две проблемы:

Доступны только биномиальная и Пуассона (проблема поменьше).
В отличие от остального statsmodels, они основаны на байесовской статистике (проблема побольше).

Байесовская статистика – это совершенно особый зверь как с позиции своих алгоритмов, так и с философской точки зрения. Более-менее нормальное освещение этой темы – отдельный разговор, поэтому сейчас мы не будем разбираться, как делать GLMM.

Хорошая новость: в GEE этот функционал реализован так же, как и в остальном statsmodels! Так что теперь у вас есть дополнительная, глубоко прагматическая причина пользоваться GEE в подобных ситуациях. Вдобавок к этому есть отдельные классы для мультиномиальной регрессии и даже порядковой регрессии! Заставляет задуматься, нет ли у авторов пакета персонального фаворитизма в сторону GEE.

Построение таких моделей практически идентично классическим GLM: мы создаём класс с нужным нам распределением отклика и подаём его в качестве аргумента family. Попробуем, например, оценить связь социально-экономического статуса с тем, является ли человек нацменьшинством (minority). Это классическая проблема логистической регрессии и реализуется она следующим кодом:

fam = sm.families.Binomial() #класс с биномиальным распределением

gee_log = smf.gee('minority ~ ses',
        data = classes, groups = classes['schoolid'],
        cov_struct = ind, family = fam).fit()
        
gee_log.summary()

Как видно из таблицы, взаимосвязь отрицательная – более низкий социальный статус скорее характерен для нацменьшинств:

Альтернатива statsmodels для LMM

Всё это время мы использовали statsmodels для построения моделей со смешанными эффектами, однако в Python существует более специализированный пакет под названием pymer4. По ряду причин он превосходит statsmodels:

Более удобный и гибкий синтаксис построения моделей с помощью тех же формул
Поддержка GLM (распределений не очень много, но они самые распространённые)
Возможность оформлять результаты моделей как дисперсионный анализ со всеми сопутствующими фичами
Ряд вспомогательных статистических процедур для оценки значимости результатов
Возможность рисовать получившиеся модели нативными функциями

Иными словами, для (G)LMM pymer4 является крайне удобным и гибким инструментом. Если вы заставите его работать.

Дело в том, что pymer4 полагается на код из языка R и требует для этого установку пакета rpy2. Последний печально известен тем, что заставить его работать не всегда просто, а если у вас Windows – очень сложно, если вообще возможно. По этой причине мы и не используем здесь этот пакет. Вряд ли этот туториал будет вам полезен, если вы не сможете пользоваться этим кодом, верно же? Тем не менее, если вам всё же выйдет с ним справиться, то очень рекомендуем опробовать его на практике. Ну и, конечно, ознакомиться с неповторимым оригиналом.