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

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

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

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

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

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

Multilevel models are regressions with coefficients that can vary by groups
Gelman, A. and Hill, J. (2006). Data Analysis Using Regression and Multilevel/Hierarchical Models [1]

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

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

Парадокс Симпсона. Автор: Schutz - собственная работа, Общественное достояние, Ссылка
Привидите собственный пример данных, имеющих многоуровневую структуру.

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

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

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

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

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

$$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). Основной предиктор уровня радона - это этаж, на котором проводилось измерение: подвал или первый этаж. Поскольку радон содержится в грунте, чем ближе к земле проиходит измерение, тем выше должен быть его уровень. Этот предиктор был измерен на уровне домохозяйства, т.е. на первом уровне. На уровне окруда есть дургой важных предиктор - средний уровень радона в данной местности. Задачей исследования, таким образом, будет отследить влияние предикторов на перовом и втором уровне на содержание радона в здании.

radon = read.csv("https://raw.githubusercontent.com/pymc-devs/pymc3/master/pymc3/examples/data/radon.csv")
radon$county = as.factor(radon$county)
dim(radon)
  1. 919
  2. 30
head(radon)
Xidnumstatestate2stfipszipregiontypebldgfloorroompcterradjwtdupflagzipflagcntyfipscountyfipsUppmcounty_codelog_radon
0 5081 MN MN 27 55735 5 1 1 3 9.7 1146.4992 1 0 1 AITKIN 27001 0.502054 0 0.83290912
1 5082 MN MN 27 55748 5 1 0 4 14.5 471.3662 0 0 1 AITKIN 27001 0.502054 0 0.83290912
2 5083 MN MN 27 55748 5 1 0 4 9.6 433.3167 0 0 1 AITKIN 27001 0.502054 0 1.09861229
3 5084 MN MN 27 56469 5 1 0 4 24.3 461.6237 0 0 1 AITKIN 27001 0.502054 0 0.09531018
4 5085 MN MN 27 55011 3 1 0 4 13.8 433.3167 0 0 3 ANOKA 27003 0.428565 1 1.16315081
5 5086 MN MN 27 55014 3 1 0 4 12.8 471.3662 0 0 3 ANOKA 27003 0.428565 1 0.95551145

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

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

library("lme4")

## Varying intercept & slopes w/ no group level predictors
M0 <- lmer(log_radon ~ floor + (1 | county), data=radon)

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

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

summary(M0)
Linear mixed model fit by REML ['lmerMod']
Formula: log_radon ~ floor + (1 | county)
Data: radon

REML criterion at convergence: 2097.5

Scaled residuals:
Min 1Q Median 3Q Max
-4.3597 -0.6253 -0.0016 0.6389 3.4965

Random effects:
Groups Name Variance Std.Dev.
county (Intercept) 0.09948 0.3154
Residual 0.52676 0.7258
Number of obs: 919, groups: county, 85

Fixed effects:
Estimate Std. Error t value
(Intercept) 1.49239 0.04955 30.117
floor -0.66289 0.06765 -9.798

Correlation of Fixed Effects:
(Intr)
floor -0.288

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

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

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

head(coef(M0)$county)  # для первых шести групп из переменной группового уровня county
(Intercept)floor
AITKIN1.2291388 -0.6628887
ANOKA0.9811502 -0.6628887
BECKER1.5066938 -0.6628887
BELTRAMI1.5377713 -0.6628887
BENTON1.4732718 -0.6628887
BIG STONE1.5084495 -0.6628887

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

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

head(ranef(M0)$county)
(Intercept)
AITKIN-0.26325439
ANOKA-0.51124297
BECKER 0.01430059
BELTRAMI 0.04537812
BENTON-0.01912135
BIG STONE 0.01605636

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

fixef(M0)
(Intercept)
1.49239317405761
floor
-0.662888725670425
(ranef(M0)$county[, 1]  + fixef(M0)['(Intercept)'][[1]])[1:5]
  1. 1.22913878203744
  2. 0.981150201857783
  3. 1.50669376251355
  4. 1.53777129831703
  5. 1.47327182833173
# вспомогательный код для отрисовки всех графиков в SVG
library("Cairo")
options(jupyter.plot_mimetypes = 'image/svg+xml')
options(warn=-1) # ALERT: отключаем предупреждения из-за кириллицы глобально

# предсказываем уровень радона и визуализируем
pr0 <- predict(M0, radon)
ggplot(data=radon, aes(floor, pr0)) + geom_line(aes(color = county), show.legend = FALSE) +
    labs(x="Этаж (подвал или 1-й этаж)", y="log(Радон)")

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

library("sjPlot")
set_theme(axis.textsize = .7)
plot_model(M0, type = "re")

Линиями на этом графиками обозначины доверительные интервалы для величины случайного эффекта в группах. Если какая-то линия не пересекает вертикальную прямую, соответствующую 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].

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

var_f <- var(predict(M0, re.form=NA))

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

(vars <- as.data.frame(VarCorr(M0)))
grpvar1var2vcovsdcor
county (Intercept)NA 0.09948229 0.3154081
Residual NA NA 0.52675634 0.7257798
var_l <- vars1[1,4]
var_e <- vars1[2,4]
varf/(varf+var_l+var_e)
0.0888194098019492
(var_f+var_l)/(var_f+var_l+var_e)
0.233566679676481

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

tab_model(M0)
Computing p-values via Wald-statistics approximation (treating t as Wald z).
  log radon
Predictors Estimates CI p
(Intercept) 1.49 1.40 – 1.59 <0.001
floor -0.66 -0.80 – -0.53 <0.001
Random Effects
σ2 0.53
τ00 county 0.10
ICC county 0.16
Observations 919
Marginal R2 / Conditional R2 0.089 / 0.234

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

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

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

M1 <- lmer(log_radon ~ floor + (floor - 1 | county), data=radon)
summary(M1)
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: log_radon ~ floor + (floor - 1 | county)
Data: radon

AIC BIC logLik deviance df.resid
2177.0 2196.2 -1084.5 2169.0 915

Scaled residuals:
Min 1Q Median 3Q Max
-3.9444 -0.6790 0.0306 0.6896 3.2931

Random effects:
Groups Name Variance Std.Dev.
county floor 0.0924 0.3040
Residual 0.6081 0.7798
Number of obs: 919, groups: county, 85

Fixed effects:
Estimate Std. Error t value
(Intercept) 1.36241 0.02818 48.355
floor -0.53455 0.08396 -6.367

Correlation of Fixed Effects:
(Intr)
floor -0.336
head(coef(M1)$county)
(Intercept)floor
AITKIN1.36241 -0.5338824
ANOKA1.36241 -0.8334118
BECKER1.36241 -0.5190840
BELTRAMI1.36241 -0.5458936
BENTON1.36241 -0.4988334
BIG STONE1.36241 -0.5345482
predict(M1, radon)[1]
1: 0.828527444546068
pr1 <- predict(M1, radon)
ggplot(data=radon, aes(floor, pr1)) + geom_line(aes(color = county), show.legend = FALSE) +
    labs(x="Этаж (подвал или 1-й этаж)", y="log(Радон)")
Модель с неизменным свободным коэффициентом, но изменяющимся углом наклон довольно редко используется, но иногда она незаменима. Приведите пример такой ситуации.

Модель со случайным свободным коэффициентом и углом наклона (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}.$$

M2 <- lmer(log_radon ~ floor + (floor|county), data=radon)
tab_model(M2)
Computing p-values via Wald-statistics approximation (treating t as Wald z).
Caution! ICC for random-slope-intercept models usually not meaningful. Use `adjusted = TRUE` to use the mean random effect variance to calculate the ICC. See 'Note' in `?icc`.
  log radon
Predictors Estimates CI p
(Intercept) 1.49 1.39 – 1.60 <0.001
floor -0.65 -0.82 – -0.49 <0.001
Random Effects
σ2 0.51
τ00 county 0.11
τ11 county.floor 0.11
ρ01 county -0.37
ICC county 0.18
Observations 919
Marginal R2 / Conditional R2 0.086 / 0.257
pr2 <- predict(M2, radon)
set_theme(axis.textsize = 1)
ggplot(data=radon, aes(floor, pr2)) + geom_line(aes(color = county), show.legend = FALSE) +
    labs(x="Этаж (подвал или 1-й этаж)", y="log(Радон)")
set_theme(axis.textsize = .7)
plot_model(M2, type = "re")