Дата занятий: 24 октября

1. Постановка задачи и данные

Рассмотрим пример оценивания моделей со смешанными эффектами на кросс-секционных данных из статьи [Luke D.A., Krauss M., 2004], которую Вы можете скачать по ссылке. Авторы изучают факторы голосования членов Конгресса США в интересах табачной индустрии. В фокусе внимания – различия во взаимосвязи размера финансирования комитетами политического действия (PAC) от табачных корпораций и голосования в интересах табачной индустрии. В частности, авторы тестируют, есть ли различия в этой взаимосвязи между представителями Демократической и Республиканской партий. Воспроизведем часть результатов, полученных в исследовании, а в домашнем задании Вас ожидает продолжение.

Ниже представлено краткое описание переменных:

  • state - штат (группирующая переменная)
  • lastname - член Конгресса

Переменные на первом уровне (изменяются по членам Конгресса):

  • votepct - доля голосов, отданных членом Конгресса в поддержку табачной индустрии (зависимая переменная)
  • party - дамми-переменная: представитель Демократической или Республиканской партии (1 - Республиканская партия)
  • money - размер финансирования от PAC табачных корпораций (в тысячах долларов)

Предиктор на втором уровне (изменяется по штатам, при этом внутригрупповой изменчивости нет):

  • acres - площадь табачных плантаций (в тысячах акров)

Откроем массив данных, группирующую переменную (штат) сразу для последующего анализа зададим как факторную.

ME<-read_dta("C:/Users/User/OneDrive/Рабочий стол/RAPDC_lab3.dta")

ME$state_id <- as.factor(ME$state) 
ME <- na.omit(ME)

head(ME)

Как мы видим, данные являются вложенными: более крупные единицы анализа – штаты (2-ой уровень), внутри которых наблюдается изменчивость по членам Конгресса (1-ый уровень). При этом отсылка к такой структуре данных не может быть единственным основанием для оценивания моделей со смешанными эффектами. Не забываем также о том, что для оценивания ME-моделей на втором уровне должно быть достаточно большое количество наблюдений (минимальное количество – порядка 30 для самой базовой экономной спецификации модели, с ростом количества параметров растут требования и к объему выборки). В нашем случае – 50 штатов на втором уровне. Для более полного обоснования мы отмечаем, что интерес представляют различия во взаимосвязи, а не только разные стартовые значения голосования в интересах табачной индустрии: а именно, авторы статьи хотели проверить, есть ли между штатами различия во взаимосвязи финансирования комитетами политического действия от табачных корпораций и голосования в интересах табачной индустрии. Проиллюстрируем различия между штатами на графике ниже:

ggplot(ME, aes(x=money, y=votepct)) + 
  geom_smooth(method = "lm", se = FALSE) +
  geom_point() +
  facet_wrap("state_id")

2. Эффект сжатия

Видно, что данные несбалансированные, о чем говорит разное количество наблюдений в штатах. Модель со смешанными эффектами несколько сглаживает ситуацию несбалансированности. Здесь надо вспомнить про так называемый “эффект сжатия” (shrinkage effect). С одной стороны, считается среднее общее по всей выборке (будь то константа или взаимосвязь в среднем) – отсылка к фиксированному эффекту, с другой стороны, считается среднее по j-ой подгруппе (то есть, именно специфика j-ого штата). Если по штату (в более общем смысле – пространственной единице на втором уровне) мало наблюдений, то значения (константы / коэффициента взаимосвязи) для такого штата будет “стягиваться” к среднему общему по всей выборке в силу того, что такой штат рассматривается в рамках ME-модели как “ненадежный”, иными словами, на такую информацию в условиях ограниченного количества наблюдений нельзя полагаться. Кстати, интерактивные картинки того, как происходит “сжатие”, можно посмотреть здесь.

3. Bottom-up strategy

Реализуем стратегию bottom-up при уточнении спецификации модели со смешанными эффектами, то есть, начинаем с простейшей ANOVA-модели, а далее будем последовательно усложнять модель.

1. ANOVA-модель

Оценим нулевую модель (то есть, без объясняющих переменных) и рассчитаем на основе полученной информации ICC.

Пояснение по коду: при использовании lmer для оценивания ME-модели, до скобок записываем фиксированные эффекты, к примеру, в нашем случае используется фиксированный эффект только на константу, поэтому ставим 1; далее в скобках записываем случайные эффекты, так, 1 в скобках означает случайный эффект на константу. Не забываем после вертикальной разделительной черты указать группирующую переменную – штат. Через логический оператор обозначаем, какой метод оценивания используем: REML = FALSE означает, что мы отменяем оценивание через restricted maximum likelihood method (подходит, если у Вас количество наблюдений на втором уровне на грани – 30 или чуть меньше), вместо этого предпочитаем классический MLE

null <- lmer(votepct ~ 1 + (1|state_id), REML = FALSE, data = ME)
summary(null)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ 1 + (1 | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##    319.7    332.5   -156.9    313.7      524 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23742 -0.85315  0.03887  0.75955  2.10386 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  state_id (Intercept) 0.03534  0.1880  
##  Residual             0.09264  0.3044  
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.53094    0.03114 45.57848   17.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

В выдаче мы видим две части: случайные эффекты (там вариация и для \(u_{oj}\), и для остатков на первом уровне – residuals) и фиксированные эффекты (оценка коэффициента константы, составляющая 0.53094 – среднее значение доли голосов в поддержку табачной индустрии по всей выборке). На основе информации, представленной в части Random effects, можем рассчитать ICC, поделив вариацию зависимой переменной на уровне штатов на общую вариацию.

0.03534 / (0.03534+0.09264) 
## [1] 0.2761369

Получается, что около 28% изменчивости доли голосов в поддержку табачной индустрии объясняется на уровне штатов, оставшаяся изменчивость – на уровне членов Конгресса.

Можем также получить ICC автоматически следующим образом:

performance::icc(null)
## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.276
##   Conditional ICC: 0.276

Здесь скорректированное значение ICC и без корректировки совпадают, но так бывает не всегда. В дальнейшем Вам понадобится скорректированный ICC.

Казалось бы, все просто: видим значение 0.28 в качестве ICC – делаем вывод о том, что нам нужен такой класс моделей, который работает с неоднородными данными. Однако не стоит торопиться. Вполне может быть так, что значение ICC (примерно 0.28) – результат наличия двух-трех нетипичных наблюдений. Во избежание ошибочной интерпретации и поспешного вывода мы опираемся на следующий график, сравнивающий средние наблюдаемые (те, что в наших данных) с результатами группировки на основе многократного извлечения случайных выборок из имеющейся совокупности данных. Давайте разъясним, что происходит пошагово, и как это отражается на графике ниже.

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

  2. Генерим случайным образом подгруппы (псевдо-штаты) из имеющейся совокупности данных и считаем средние по полученным подгруппам, упорядочиваем от меньшего к большему. При генерации мы сохраняем количество подгрупп (то есть, в нашем случае – 50 по количеству штатов) и размер исходных выборок.

  3. Много раз повторяем процедуру из второго шага. Если Вы посмотрите на код ниже, то сейчас мы задаем 1000 репликаций. Далее усредняем полученные результаты по тысячи репликациям и получаем сглаженную линию с доверительным интервалом, рассчитанным на основе указанной процедуры бутстрэпа (см. график).

graph.ran.mean(ME$votepct, ME$state_id, nreps=1000, bootci=TRUE)

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

Выгрузку предсказанных значений случайных эффектов (BLUP) можно получить следующим образом:

ranef(null)
## $state_id
##    (Intercept)
## AK  0.11862432
## AL  0.22093945
## AR  0.07706954
## AZ  0.11410831
## CA -0.08432885
## CO  0.11590471
## CT -0.27682984
## DE -0.08694932
## FL  0.03672017
## GA  0.15004250
## HI -0.19312165
## IA -0.01096309
## ID  0.17824836
## IL -0.03809144
## IN  0.12170401
## KS  0.13242530
## KY  0.30474305
## LA  0.14575977
## MA -0.41172037
## MD -0.16884885
## ME -0.15538418
## MI -0.11329025
## MN -0.11869996
## MO  0.09653619
## MS  0.13394838
## MT  0.01846444
## NC  0.25355886
## ND -0.20275882
## NE  0.08519685
## NH  0.14588531
## NJ -0.18059014
## NM  0.01350197
## NV -0.15383526
## NY -0.17127608
## OH -0.03952206
## OK  0.26547879
## OR -0.15167174
## PA -0.05883805
## RI -0.24919611
## SC  0.18945244
## SD -0.13601415
## TN  0.15722713
## TX  0.06704161
## UT  0.03820302
## VA  0.12892897
## VT -0.13544721
## WA -0.12374887
## WI -0.11458797
## WV -0.12873965
## WY  0.19474048
## 
## with conditional variances for "state_id"

Так, для Аляски (AK) BLUP для константы составляет 0.1186, то есть, среднее значение доли голосов в поддержку табачной индустрии для Аляски больше на 0.1186, чем среднее по всей выборке.

2. Добавляем фиксированные эффекты на предикторы

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

options(scipen = 999)
model1 <- lmer(votepct ~ money + (1|state_id), REML = FALSE, data = ME)
summary(model1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + (1 | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##    237.6    254.7   -114.8    229.6      523 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.33857 -0.82138 -0.07304  0.78857  2.21384 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  state_id (Intercept) 0.01707  0.1307  
##  Residual             0.08208  0.2865  
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##                Estimate  Std. Error          df t value            Pr(>|t|)    
## (Intercept)   0.4287491   0.0259064  52.4900288  16.550 <0.0000000000000002 ***
## money         0.0076100   0.0007771 500.1512953   9.793 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## money -0.404

В этой модели по-прежнему только случайный эффект на константу. Давайте изобразим на графике индивидуальные стартовые точки (средние значения доли голосов в поддержку табачной индустрии при условии отсутствия финансирования) для каждого штата, что представляет сумму фиксированного эффекта на константу и случайного эффекта (\(\gamma_{00} + u_{oj}\)), также на графике будут изображены доверительные интервалы для константы.

dotplot(ranef(model1, condVar=TRUE))
## $state_id

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

В ряде случаев не очень удобно интерпретировать результаты в терминах отклонений, вместо этого хотелось бы сразу получить среднее стартовое значение и взаимосвязь ключевого предиктора и зависимой переменной в j-ой пространственной единице. Рассчитаем оценки констант и коэффициента взаимосвязи финансирования от PAC табачных корпораций и доли голосов в поддержку табачной индустрии для каждого штата. В квадратных скобках обозначим номер параметра, который нас интересует. [,1] – для константы, [,2] – для коэффициента взаимосвязи.

model1_intercept <- coef(model1)$state_id[,1] 
model1_slope <- coef(model1)$state_id[,2] 
state_unique <- unique(ME$state_id)
data.frame(state_unique, model1_intercept, model1_slope)

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

Продемонстрируем, как можно сравнивать модели. На данном этапе сравним нулевую модель и первую модель. Данные модели являются вложенными: можем получить model1 из нулевой модели посредством добавления параметра – фиксированного эффекта для money. Проверяем гипотезу о равенстве коэффициента при money нулю. Или иными словами, в рамках теста отношения правдоподобия для вложенных моделей мы задаем вопрос, а вносит ли предиктор (money) значимый вклад в объяснение изменчивости доли голосов в поддержку табачной индустрии?

anova(null, model1)

В выдаче смотрим на статистику \(\chi^2\), p-value близко к 0, следовательно, нулевую гипотезу отвергаем. То есть, выбираем модель менее экономную, с включенным предиктором money. Также можно сравнить информационные критерии, показывающие, сколько информации мы теряем. Важно найти баланс между нагруженностью модели и ее информативностью. BIC больше штрафует за добавление параметров по сравнению с AIC. Как по AIC, так и по BIC мы видим, что лучше model1 (так как меньше значения AIC, BIC).

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

model2 <- lmer(votepct ~ money + party + (1|state_id), REML = FALSE, data = ME)
summary(model2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + party + (1 | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##   -302.3   -280.9    156.1   -312.3      522 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3363 -0.5885 -0.0693  0.5980  4.1462 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  state_id (Intercept) 0.006192 0.07869 
##  Residual             0.029328 0.17125 
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##                Estimate  Std. Error          df t value            Pr(>|t|)    
## (Intercept)   0.1997207   0.0172421  76.4322756  11.583 <0.0000000000000002 ***
## money         0.0044760   0.0004758 507.4162925   9.407 <0.0000000000000002 ***
## party         0.4961530   0.0161260 519.2072838  30.767 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) money 
## money -0.262       
## party -0.432 -0.214
anova(model1, model2)

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

На этом не останавливаемся и добавляем фиксированный эффект для переменной площадь табачных плантаций – характеристике на втором уровне (есть только межгрупповая изменчивость). Здесь мы следуем логике авторов, указывающих на необходимость проконтролировать переменные – характеристики штатов, а не только характеристики членов Конгресса. Мы видим, что у оценки коэффициента при acres пограничный уровень значимости (тест отношения правдоподобия, сравнивающий модель 2 и модель 3, показывает схожие результаты), однако из этого пока не следует, что этот предиктор нужно сразу удалять из модели. Вспомним, что у авторов статьи также была идея протестировать, есть ли различия во взаимосвязи характеристик членов Конгресса и доли голосов в поддержку табачной индустрии в зависимости от площади табачных плантаций, вполне может быть так, что acres окажется значимой переменной-модератором.

model3 <- lmer(votepct ~ money + party + acres + (1|state_id), REML = FALSE, data = ME)
summary(model3)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + party + acres + (1 | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##   -304.1   -278.5    158.1   -316.1      521 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3490 -0.5823 -0.0572  0.5970  4.1815 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  state_id (Intercept) 0.005664 0.07526 
##  Residual             0.029258 0.17105 
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##                Estimate  Std. Error          df t value             Pr(>|t|)
## (Intercept)   0.1929140   0.0171814  72.2590879  11.228 < 0.0000000000000002
## money         0.0041863   0.0005002 512.9078659   8.369 0.000000000000000554
## party         0.4970312   0.0160821 520.4813746  30.906 < 0.0000000000000002
## acres         0.0006535   0.0003293  49.2721164   1.985               0.0528
##                
## (Intercept) ***
## money       ***
## party       ***
## acres       .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) money  party 
## money -0.187              
## party -0.434 -0.207       
## acres -0.186 -0.323  0.014
anova(model2, model3)

3. Добавляем случайные эффекты для того, чтобы учесть различия во взаимосвязи

На предыдущем шаге мы не учли возможные различия между штатами во взаимосвязи характеристик членов Конгресса и доли голосов в поддержку табачной индустрии. Исправим это упущение в последующих спецификациях. Так, в модели 4.1 учитывается различия во взаимосвязи финансирования от PAC табачных корпораций и отклика – для этого добавляем \(u_{1j}\).

model4.1 <- lmer(votepct ~ money + party + acres + (1 + money|state_id), REML = FALSE, data = ME)
summary(model4.1) # non-convergence
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + party + acres + (1 + money | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##   -305.0   -270.9    160.5   -321.0      519 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2680 -0.6009 -0.0617  0.5485  4.0688 
## 
## Random effects:
##  Groups   Name        Variance    Std.Dev. Corr 
##  state_id (Intercept) 0.008189956 0.090498      
##           money       0.000003755 0.001938 -0.71
##  Residual             0.028305100 0.168241      
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##                Estimate  Std. Error          df t value             Pr(>|t|)
## (Intercept)   0.1982417   0.0187228  48.4372247  10.588   0.0000000000000337
## money         0.0044453   0.0006069  31.6784598   7.325   0.0000000266203122
## party         0.4927196   0.0161233 519.0842813  30.559 < 0.0000000000000002
## acres         0.0007538   0.0003089   9.7119050   2.440               0.0355
##                
## (Intercept) ***
## money       ***
## party       ***
## acres       *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) money  party 
## money -0.392              
## party -0.413 -0.179       
## acres -0.131 -0.300  0.022
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0597448 (tol = 0.002, component 1)

Что-то пошло не так, а именно, мы видим предупреждение о несходимости алгоритма.

model4.1 <- lmer(votepct ~ money + party + acres + (1 + money|state_id), REML = FALSE, data = ME)
summary(model4.1) # non-convergence
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + party + acres + (1 + money | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##   -305.0   -270.9    160.5   -321.0      519 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2680 -0.6009 -0.0617  0.5485  4.0688 
## 
## Random effects:
##  Groups   Name        Variance    Std.Dev. Corr 
##  state_id (Intercept) 0.008189956 0.090498      
##           money       0.000003755 0.001938 -0.71
##  Residual             0.028305100 0.168241      
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##                Estimate  Std. Error          df t value             Pr(>|t|)
## (Intercept)   0.1982417   0.0187228  48.4372247  10.588   0.0000000000000337
## money         0.0044453   0.0006069  31.6784598   7.325   0.0000000266203122
## party         0.4927196   0.0161233 519.0842813  30.559 < 0.0000000000000002
## acres         0.0007538   0.0003089   9.7119050   2.440               0.0355
##                
## (Intercept) ***
## money       ***
## party       ***
## acres       *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) money  party 
## money -0.392              
## party -0.413 -0.179       
## acres -0.131 -0.300  0.022
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0597448 (tol = 0.002, component 1)
model4.1.1 <- lmer(votepct ~ money + (1 + money|state_id), REML = FALSE, data = ME)
## boundary (singular) fit: see help('isSingular')
summary(model4.1.1) # non-convergence
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + (1 + money | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##    223.9    249.5   -105.9    211.9      521 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.05956 -0.76422 -0.08976  0.73254  2.35112 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev. Corr 
##  state_id (Intercept) 0.03470289 0.186287      
##           money       0.00001856 0.004308 -1.00
##  Residual             0.07775704 0.278849      
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##               Estimate Std. Error         df t value             Pr(>|t|)    
## (Intercept)  0.4465684  0.0320924 41.5101980  13.915 < 0.0000000000000002 ***
## money        0.0081311  0.0009472 30.8182293   8.585        0.00000000113 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## money -0.791
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

В чем может быть проблема? Вариация случайного эффекта для money ничтожна мала. Давайте посмотрим на следующую картинку: сравним предсказание в среднем по всей выборке и предсказание с учетом специфики по штатам (случайного эффекта как на константу, так и на наклон). Команда predict(model) по умолчанию рассчитывает предсказанные значения с учетом специфики по подгруппам. Если используем predict(model, re.form = NA), то предсказание создается только на основе фиксированных эффектов, и как следствие, одинаково для всех штатов.

ME %>% 
  # add predicted values to the dataset
  mutate(pr_re = predict(model4.1.1), pr_fe = predict(model4.1.1, re.form = NA)) %>% 
  # visualize
  ggplot(aes(x=money, y=pr_re, group = state_id)) + theme_light() + geom_line() +
  geom_line(color = "red", aes(money, pr_fe)) + facet_wrap(~state_id)

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

Что если по штатам различается эффект партийной принадлежности? Проверим аналогичной спецификацией, добавив случайный эффект на party.

model4.2 <- lmer(votepct ~ party + (1 + party|state_id), REML = FALSE, data = ME)
summary(model4.2) 
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ party + (1 + party | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##   -242.4   -216.8    127.2   -254.4      521 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7731 -0.5559  0.0134  0.5491  4.7581 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  state_id (Intercept) 0.02714  0.1647        
##           party       0.01063  0.1031   -0.88
##  Residual             0.03021  0.1738        
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  0.25694    0.02800 40.00208   9.177      0.0000000000216 ***
## party        0.51259    0.02386 36.35950  21.483 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## party -0.790

И модель с полным набором исходных предикторов:

model4.2.1 <- lmer(votepct ~ money + acres + party + (1 + party|state_id), REML = FALSE, data = ME)
summary(model4.2.1) 
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + acres + party + (1 + party | state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##   -312.0   -277.8    164.0   -328.0      519 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7460 -0.5764 -0.0002  0.5632  4.5800 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  state_id (Intercept) 0.010816 0.10400       
##           party       0.008485 0.09211  -0.73
##  Residual             0.027284 0.16518       
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##                Estimate  Std. Error          df t value             Pr(>|t|)
## (Intercept)   0.2047760   0.0209332  34.7334789   9.782  0.00000000001630022
## money         0.0040650   0.0004947 510.9469193   8.217  0.00000000000000173
## acres         0.0005924   0.0003174  36.9971747   1.866               0.0699
## party         0.4882487   0.0222164  37.8990776  21.977 < 0.0000000000000002
##                
## (Intercept) ***
## money       ***
## acres       .  
## party       ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) money  acres 
## money -0.164              
## acres -0.142 -0.325       
## party -0.669 -0.137 -0.001
dotplot(ranef(model4.2.1, condVar=TRUE))
## $state_id

В выдаче модели 4.2.1 в части случайных эффектов есть оценка такого параметра Corr, составляющая \(-0.73\). Это корреляция между случайными эффектами для константы и взаимосвязи партийной принадлежности и доли голосов в поддержку табачной индустрии. Изобразим на графике эту взаимосвязь. Получается, что чем ниже исходный уровень голосования в поддержку табачной индустрии, тем менее выражена взаимосвязь партийной принадлежности и данной зависимой переменной. С одной стороны, за этим может стоять содержательная история, которую интересно “раскрутить” в исследовании. С другой стороны, значимую корреляцию на данном этапе мы можем получить и в результате некоторых пропущенных переменных.

random4.2.1 <- ranef(model4.2.1)
unique <- cbind((random4.2.1$state_id[,1]),(random4.2.1$state_id[,2]))
plot(unique, 
     xlab = 'Intercepts', ylab = "Slopes ",
     main = "Correlation between BLUP Intercepts and Slopes")

Представим, что корреляция не так уж велика (к примеру, по модулю не превышает 0.3). Тогда можно пожертвовать этим параметром, и во имя более экономной спецификации внести ограничение, заключающееся в том, что данная корреляция равна 0. Такое ограничение можно прописать в явном виде, изменив структуру корреляционной матрицы. Из library multilevel используем lme для переоценки модели, при этом предлагаем диагональную матрицу pdDiag для случайных эффектов: в такой матрице элементы вне главной диагонали равны 0.

model_nlme_1 <- lme(votepct ~  money + party + acres, random = list(state_id = pdDiag(~1 + ~party)), method = "ML", data = ME)
summary(model_nlme_1)
## Linear mixed-effects model fit by maximum likelihood
##   Data: ME 
##         AIC       BIC   logLik
##   -307.9102 -278.0398 160.9551
## 
## Random effects:
##  Formula: ~1 + ~party | state_id
##  Structure: Diagonal
##         (Intercept)      party  Residual
## StdDev:  0.07850753 0.06907026 0.1661866
## 
## Fixed effects:  votepct ~ money + party + acres 
##                 Value   Std.Error  DF   t-value p-value
## (Intercept) 0.1934876 0.017789906 475 10.876257  0.0000
## money       0.0041458 0.000497678 475  8.330253  0.0000
## party       0.4858930 0.019899316 475 24.417574  0.0000
## acres       0.0009781 0.000365842  48  2.673609  0.0102
##  Correlation: 
##       (Intr) money  party 
## money -0.185              
## party -0.382 -0.153       
## acres -0.204 -0.298 -0.006
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.706280444 -0.602561498 -0.003090901  0.538173214  4.481036555 
## 
## Number of Observations: 527
## Number of Groups: 50

Несмотря на то, что модель со смешанными эффектами в явном виде учитывает специфику единиц анализа на втором уровне (в нашем случае – штатов), это не избавляет нас от необходимости выявления влиятельных наблюдений. По графику ниже на основании меры Кука видно, что такими влиятельными наблюдениями для модели 4.2 являются штаты Кентукки, Нью-Джерси и Калифорния.

inf <- influence(model4.2, group = "state_id", data = ME)
cooks.distance.estex(inf, sort=TRUE)
##             [,1]
## AK 0.00004054182
## NH 0.00017916274
## AZ 0.00075076707
## HI 0.00084344245
## MN 0.00104960260
## WV 0.00187762982
## KS 0.00221803193
## NE 0.00263385176
## MT 0.00272537048
## ND 0.00344670156
## PA 0.00418566106
## FL 0.00453735574
## AR 0.00455857983
## NV 0.00456286730
## ID 0.00518232184
## OH 0.00519854919
## OR 0.00616596050
## WY 0.00791596535
## MI 0.00839486907
## SD 0.00841956199
## IA 0.00889327204
## WA 0.01177476240
## NM 0.01235928103
## TX 0.01254590812
## OK 0.01283935538
## MS 0.01516498980
## CO 0.02028902789
## IN 0.02223470186
## MA 0.02245247585
## LA 0.02251916324
## MD 0.02253094147
## VT 0.02463733979
## MO 0.02611022862
## RI 0.02679411108
## NY 0.02710795542
## SC 0.02892248805
## VA 0.02923190794
## DE 0.03121342084
## GA 0.03553045880
## WI 0.03696959394
## UT 0.03831808836
## IL 0.04387885446
## CT 0.05006207459
## ME 0.06476927772
## TN 0.07166585615
## AL 0.07237961159
## NC 0.07920816415
## KY 0.08540665468
## CA 0.09817495570
## NJ 0.30782079613
plot(inf, which = "cook", xlab = "COOK'S MEASURE", cutoff = 4/length(state_unique))

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

plot_model(model4.2.1, type = "diag")[1]
## [[1]]

plot_model(model4.2.1, type = "diag")[3]
## [[1]]

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

plot_model(model4.2.1, type = "diag")[2]
## [[1]]
## [[1]]$state_id

4. Добавляем cross-level interaction

В соответствии с одной из гипотез авторов добавляем в модель переменную взаимодействия между характеристиками разных уровней – партийной принадлежности и площади табачных плантаций. Оценка коэффициента при данной переменной взаимодействия оказывается значимой, что говорит о том, что в штатах с большей площадью табачных плантаций эффект партийной принадлежности на долю голосов в поддержку табачной индустрии становится менее выраженным. Оценка коэффициента при acres теперь приобретает значимость, что неудивительно, так как теперь у этого коэффициента меняется интерпретация: в данном случае интерпретируем в терминах условного эффекта, а именно, эффекта для подвыборки представителей Демократической партии (party = 0).

model5 <- lmer(votepct ~ money + party + acres + party*acres + (1 + party|state_id), REML = FALSE, data = ME)
summary(model5)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: votepct ~ money + party + acres + party * acres + (1 + party |  
##     state_id)
##    Data: ME
## 
##      AIC      BIC   logLik deviance df.resid 
##   -322.8   -284.4    170.4   -340.8      518 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7296 -0.5834  0.0081  0.5848  4.5818 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  state_id (Intercept) 0.009137 0.09559       
##           party       0.004629 0.06804  -0.71
##  Residual             0.027154 0.16479       
## Number of obs: 527, groups:  state_id, 50
## 
## Fixed effects:
##                Estimate  Std. Error          df t value             Pr(>|t|)
## (Intercept)   0.1859884   0.0202907  38.4145889   9.166  0.00000000003272077
## money         0.0040194   0.0004896 512.5065956   8.210  0.00000000000000181
## party         0.5164640   0.0209246  38.8941096  24.682 < 0.0000000000000002
## acres         0.0019306   0.0004851  60.0507199   3.980             0.000188
## party:acres  -0.0017391   0.0004644  55.5266832  -3.745             0.000431
##                
## (Intercept) ***
## money       ***
## party       ***
## acres       ***
## party:acres ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) money  party  acres 
## money       -0.156                     
## party       -0.662 -0.159              
## acres       -0.272 -0.238  0.262       
## party:acres  0.235  0.037 -0.336 -0.770

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