Многоуровневые регрессионные модели

Краткая характеристика многоуровневых регрессионных моделей

Сегодня мы затронем крайне интересную тему смешанных (их ещё называют многоуровневыми или иерархическими) регрессионных моделей. Интересная она по нескольким причинам:

• многоуровневые регрессии — это очень мощный инструмент для проведения конфирматорных исследований и для исследований, измерения в которых производятся в каких-либо группах. С их помощью можно проверить достаточно сложные гипотезы, к которым не удасться подобрать ключик при помощи других методов;
• о них мало кто знает в современном мире курсеророждённых датасаентистов, поскольку, как я уже говорил, смешанные модели модели больше всего подходят для конфирматорынх исследований. Этот тип исследований достаточно редко встречается в индустрии, где основной задачей является предсказание целевой переменной, и усилия, соответственно, направлены на минимизацию ошибки этого предсказания;
• за ними лежит довольно простая, но красивая идея, о которой мы сейчас поговорим.

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

Другое называние этой группы моделей — многоуровневые регрессионные модели — присходит из-за того, что упомянутые ранее случайные эффекты (или коэффициенты) могут изменяться по группам.

Multilevel models are regressions with coefficients that can vary by groups

Чтобы лучше понять, что из себя представяют такие уровни или группы, представим ситуацию, когда мы исследуем ответы жителей разных стран на одни и те же вопросы. На этом премере видно, что в данных выделяется по меньшей мере два уровня наблюдений — индвидуальный и страновой, причём каждом уровне может быть свой набор переменных.

В процессе исследования мы строим регрессионную модель, в которой моделируем связь между доходом респондента и его удовлетворённостью жизнью. Можно взять и построить регрессионную модель на всех наблюдениях стразу, но, как вы, наверное, знаете, зависимость между счастьем и доходом далеко не одинакова для всех стран — латиноамериканские страны регулярно занимают достаточно высокие места в списке самых счастливых стран, хотя по рейтингу ВВП на душу неселения они занимают аналогичные места, но с конца (для того, чтобы узнать больше об этом кажущемся парадоксе, можете послушать подкаст Freakonomics). По этой причине более правильным решением было бы простроить отдельную регрессионную модель для каждой страны. Если этого не сделать, возникает риск столкнуться с парадоксом Симпсона, о котором мы говорили в предудущей части курса. Напомню, что парадокс Симпосна — это эффект, когда при наличии двух групп данных, в каждой из которых наблюдается одинаково направленная зависимость, при объединении этих групп направление зависимости меняется на противоположное, как на картинке:

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

Однако анализировать множество отдельных моделей крайне неудобно. Частичным решением можно считать добавление фиктивный переменных для каждой страны — обычно именно так и поступают.

Но что, если нам захочется проследить, пример, какие именно характеристики страны влияют связь между доходами и счастьем? Тут на помощь приходят многоуровневые регрессионные модели. В общем виде многоуровневую модель можно представить как регрессию (линейную или обобщенную линейную модель), в которой параметры - коэффициенты регрессии - сами являются вероятностной моделью. Эта модель второго уровня имеет собственные параметры - гиперпараметры модели, которые также оцениваются по данным. Именно эта особенность — моделирование различий между группами — отличает многоуровневые модели от классической регрессии.

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

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

$$y_i = \alpha + \beta x_i +\epsilon_i,$$ где в нашем случае $y_i$ — это доход $i$-ого человека, $\alpha$ и $\beta$ — коэффициенты модели, $x_i$ — уровень его религиозности, $\epsilon_i$ — ошибка предсказания.

Предположим, мы хотим заложить в модель возможность изменения среднего уровня религиозности для каждой страны $j\in[1,J]$. Для этого необходимо, чтобы свободный член уравнения, отвечающий за положение прямой относительно оси $y$, принимал разные значения для каждой из стран. переписать свободный член уравнения следующим образом:

$$y_{ij} = \alpha_{j} + \beta x_i +\epsilon_i,$$ где $$\alpha_j=a+bu_j+\eta_j.$$

То есть $$y_{ij} = a+bu_j+\eta_j + \beta x_i +\epsilon_i,$$ где $x_i$ и $u_i$ — предикторы на индивидуальном и страновом уровне соответственно, а $\epsilon_i$ и $\eta_j$ — ошибка на этих уровнях.

В более сложном случае мы предполагаем, что не только средний уровнь религиозности может различаться по странам также варьируется характер его зависимости от дохода. Формально это выражается в разном угле наклона регрессионной прямой для каждой из стран, а ещё более формально в следующих уравнениях: $$y_{ij} = \alpha_j+\beta_j x_i+\epsilon_i,$$ где $$\alpha_j=a_0+b_0u_j+\eta_{j1},$$ $$\beta_j=a_1+b_1u_j+\eta_{j2}.$$

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

Говоря о параметрах модели, следует упомнянуть о некоторых ограничениях, связанных с количеством групп второго уровня (в прошлом примере это были страны) — в лучшем случае их должно быть не менее 50 и размер каждой группы не менее 5 наблюдений. Количество групп от 26 до 49 также допустимо, но в этом случае следует валидировать модель с использованием байесовского вывода. Модель с количеством групп на втором уровне менее 26 стоит избегать.

В конце первой вводной части я сошлюсь на выступление Анатолия Карпова о применении смешанных регрессионных моделей при анализе психологических данных.

Чем можно заменить многоуровневые регрессионные модели?

Зачем всё-таки нужны многоуровневые регрессионные модели?

• Позволяют оценивать коэффициенты для каждой из групп, причём даже при малом размере группы (>5 наблюдений).
• Позволяют оценить воздействие переменной группового уровня на переменные индивидуального уровня.
• Лучше моделируют эффекты индивидуального уровня и позволяют удобно предсказывать их значения для разных групп.
• Позволяют контролировать воздействие одной переменной на другую. Одной из основных целей регрессионного анализа является оценка того, как изменяется y, когда изменяется некоторый x, а все остальные входные данные остаются постоянными. Во многих приложениях интерес представляет не общий эффект x, а то, как этот эффект варьируется по группам. В классической статистике мы можем изучить эту вариацию, используя взаимодействия: например, определенная образовательная инновация может быть более эффективной для девочек, чем для мальчиков, или более эффективной для учащихся, которые проявили больший интерес к школе перед поступлением. Многоуровневые модели также позволяют нам изучать эффекты, которые варьируются в зависимости от группы. В классической регрессии оценки различных эффектов могут быть шумными, особенно когда в группе мало наблюдений. Многоуровневое моделирование позволяет нам более точно и надёжно оценивать эти взаимодействия в условиях ограниченных данных. и т.д.

Практика

Давате передохнём от теории и перейдём к практике. Мы будем тренировать навыки построения многоуровневых регрессионных моделей на известном наборе данных об уровне радона в жилых домах. Радон образуется при распаде радия, который находится в грунте, почве и в минеральных составляющих строительных материалов. При плохом воздухообмене в квартире он имеет свойство накапливаться, и при попадании в организм вместе с вдыхаемым воздухом радон распадается на дочерние продукты, которые облучают ткани лёгких. По данным ВОЗ, воздействие этого газа является второй после курения причиной возникновения рака лёгких, этому ВОЗ основала всемирный проект, направленный на уменьшение риска заражения радоном.

Используемый нами набор данных был собран Гельманом и использовался им для иллюстрации работы многоуровневых моделей в работе, на которую я уже ссылался [1] (кстати, подписывайтесь на его блог). На этом наборе данных удобно проверять влияние различных факторов на уровень радона.

Практически все распространённые пакеты для работы с многоуровневыми моделями созданы для R. Мы будем использовать пакет lme4. Вы также можете рассмотреть другие варианты — nlme или plm. При работе в Python можно использовать имплементацию многоуровневых регрессионных моделей из statsmodels, но там имеется довольно ограниченный набор моделей и инструментов для их тестирования и визуализации.

Есть ещё stan, который тоже может моделировать многоуровневые регрессии (и, кстати, имеет интерфейс на Python), но его основная фишка не в этом, да и сам пакет заслуживает отдельной статьи. То же самое применимо к PyMC, и о нём, думаю, мы поговорим позже.

Итак, давайте взглянем на данные. Их составляют измерения радона в различных домохозяйствах. Для удобства анализа все наблюдения были сгруппированы по округам (переменная county). Основной предиктор уровня радона - это этаж, на котором проводилось измерение: подвал или первый этаж. Поскольку радон содержится в грунте, чем ближе к земле проиходит измерение, тем выше должен быть его уровень. Этот предиктор был измерен на уровне домохозяйства, т.е. на первом уровне. На уровне окруда есть дургой важных предиктор - средний уровень радона в данной местности. Задачей исследования, таким образом, будет отследить влияние предикторов на перовом и втором уровне на содержание радона в здании.

Модель со случайным свободным коэффициентом (Random intercepts model)

Построим регрессионную модель со свободным коэффициентом, которым может принимать разные значения для разных округов. Формула этой модели, напомню, выглядит следующим образом: $$y_{ij} = \alpha_{j} + \beta x_i +\epsilon_i,$$ $$\alpha_j=a+bu_j+\eta_j.$$ Для этого вызовем функцию lmer из пакета lme4, в которую запишем формулу модели.

Синтаксис формулы не сильно отличается от такового для обычных линейных моделей, но есть несколько особенностей, связанных с обозначением многоуровневости (подробности смотрите в документации, стр. 7). В примере выше часть формулы (1 | county) означает, что мы хотим сделать свободный коэффициент случайной переменной, различающейся от округа (county) к округу.

Посмотрим на вывод этой модели. Его можно получить при помощи стандартной функции summary.

Вывод модели несколько отличается от такового для обычной линейной модели. В нём есть отдельный блок для случайных и фиксированных эффектов. Раздел, посвящённый фиксированным эффектам в целом похож на вывод обычной модели. В нем находятся значения коэффициентов и информация об их значимости. В нашём случае номер этажа, на котором проводилось измерение, находится в обратной связи с уровнем радона. Поднятие на один этаж приводит из уменьшению радона на 0.66 (мне не известны единицы измерения, к тому же зависимая переменная была прологарифмирована, поэтому сложно интерпретировать это значение). Таким образом, уравнение регрессионной прямой для среднего округа записывается как $y=1.49239-0.66289\times x$.

В части вывода, посвящённой случайным эффектам, находится стандартные отклонения и дисперсия (просто стандартное отклонение в квадрате) для переменных группового уровня. Стандартное отклонение по стране, равное 0.3154 показывает, как сильно изменяется $\alpha_j$ в уравнении модели, а число 0.7258 показывает вариацию $\epsilon_j$. Сложно как-то интерпретировать это напрямую. То же самое с полем "Correlation of Fixed Effects" — его значение редко пригодится в реальном анализе.

Куда более интересно взглянуть на список случайных эффектов для каждой группы, который можно получить при помощи функций ranef или coef. Первая возвращает только значения случайных эффектов, вторая — и фиксированных тоже.

Этот вывод означает, что для округа AITKIN регрессионное уравнение выглядит как $y=1.2291388-0.6628887\times x$, для округа ANOKA как $y=0.9811502-0.6628887\times x$ и т.д.

Функция ranef показывает, на сколько выше или ниже находится регрессионная прямая в каждом округе. Например, в тех округах, где значение отрицательно — уровень радона ниже, чем в среднем.

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

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

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

В нашем случае, например, содержание радона в зданиях из города под названием St. Louis (на графике обозначен как STLOUIS, 16-й сверху) значимо меньше, чем в среднем.

Другой способ узнать значимость коэффициентов — использовать пакет lmerTest. Достаточно просто импортировать его и вывод модели, получаемый функцией summary будет расширен.

Однако сам автор пакета lme4 остерегает от того, чтобы полагаться на p-values и доверительные интералы при оценке многоуровневых регрессий, на что у него есть ряд аргументов [2].

Коэффициент детерминации $R^2$ в многоуровневых регрессионных моделях

Вы, возможно, заметили, что в выводе стандартной функции отсутствовала информация о коэффициенте детерминации $R^2$. Причина тому — сложность определения того, что называть $R^2$ для многоуровневых регрессионных моделей. Только в 2012 году Shinichi Nakagawa и Holger Schielzeth [3] предложили вариант (дополнен Paul C.D. Johnson в 2014 [4]), который стал стандартным способом решения этой проблемы. Их предложение состояло в том, чтобы рассчитывать два варианта $R^2$ — предельный и условный. Предельный $R^2$ показывает долю дисперсии, объяснённую только лишь фиксированными эффектами:

$$R^2_{marginal}=\frac{\sigma^2_f }{\sigma^2_f+\sum^{u}_{l=1}\sigma^2_l+\sigma^2_e+\sigma^2_d},$$ где $\sigma^2_f$ — дисперсия, объяснённая фиксированными коэффициентами, $u$ — количество случайных эффектов, $\sigma^2_l$ — дисперсия в группе $l$, $\sigma^2_e$ и $\sigma^2_d$ — дисперия остатков.

Условный $R^2$ показывает долю дисперсии, объяснённую как фиксированными, так и случайными эффектами:

$$R^2_{marginal}=\frac{\sigma^2_f +\sum_{l=1}^{u}\sigma^2_l}{\sigma^2_f+\sum^{u}_{l=1}\sigma^2_l+\sigma^2_e+\sigma^2_d},$$

Давайте рассчитаем предельный и условный $R^2$ самостоятельно [5].

Дисперсия, объяснённая фиксированными коэффициентами, оценивается через дисперсию прогнозов для данных, полученном без учёта случайных переменных.

Дисперсию случайных эффектов можно получить, ипользу функцию VarCorr.

Кстати говоря, в пакете sjPlot функция tab_model, которая рисует более опрятную и информативную таблицу, чем обычная summary. В ней уже показываются значения $R^2$.

Как видно, они практически совпадают с тем, что мы рассчитали вручную.

Модель со случайным углом наклона (Random Slopes Model)

$$y_{ij} = \alpha+\beta_j x_i+\epsilon_i,$$ где $$\beta_j=a_1+b_1u_j+\eta_{j2}.$$

Модель со случайным свободным коэффициентом и углом наклона (Random Intercepts and Slopes Model)

$$y_{ij} = \alpha_j+\beta_j x_i+\epsilon_i,$$ где $$\alpha_j=a_0+b_0u_j+\eta_{j1},$$ $$\beta_j=a_1+b_1u_j+\eta_{j2}.$$

Сравнение с обычной регрессионной моделью

Сравним выводы многоуровневой регрессионной модели с выводом обычной модели, в которой группирующая переменная присутствует как обычная категориальная переменная с dummy-coding, как это обычно и бывает.

Многоуровневая модель с предиктором группового уровня

Концентрация радона в воздухе зависит от геологической обстановки. Например, граниты, в которых много урана, являются активными источниками радона. В геологии измерение содержания радона в воздухе и воде даже применяется для поиска месторождений урана и тория.

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

Многоуровневая модель с эффектом взаимодействия между предиктором группового и индивидуального уровня

Проверка многоуровновости при помощи ICC

Для того, чтобы узнать, следует ли в вашем случае строит многоуровневую модель, или можно обойтись обычно, следует рассчитать коэффициент внутриклассовой корреляции (ICC). ICC вычисляется как отношение межгрупповой дисперсии к общей дисперсии (то есть к сумме межгрупповой и внутригрупповой диспесий), а его формула выглядит следующим образом: $$ICC=\frac{st.dev(intercept)^2}{st.dev(intercept)^2 + st.dev(residual)^2}=$$ $$\frac{variance (intercept)}{(variance (intercept)+variance (residual)}$$ ICC можно интерпретировать как долю дисперсии, объясняемой многоуровневостью модели, это доля групповой дисперсии в общей дисперсии. Таким образом, ICC описывает, насколько различаются группы между собой группы. ICC больше 0.1 обычно считается достаточным для построения многоуровневых моделей, хотя чётких критериев тут нет.

Обычно ICC рассчитывается для нулевой модели (её ещё назвают пустой моделью), т.е. модели, в которой заложена только группирующая переменная. Фактически, это обычная модель ANOVA.

Рассчитать ICC также можно при помощи функции icc из пакета sjstats, который помимо неё содержит множество разнообразных функций для работы с регрессионными моделями.

Как видно, эти значения почти плоностью повторяют рассчитанные вручную.

Ограничения многоуровневой регрессионной модели

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