Модерация или очищение эффекта в R на примере интроверсии замужних женщин

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

1. X - некий фактор, который имеет корреляцию с Y, и если мы создадим уравнение Y ~ X, то связь будет значима.
2. M - фактор, который коррелирует с X и Y. И если мы запишем Y ~ M и X ~ M, то в обоих случая связь будет значима. 
3. Но если мы запишем Y ~ X + M, то значимость X теряется, потому что на самом деле связь есть только у M с Y

Это и есть эффект модерации. Думаю, не нужно сильно погружаться в объяснение важности этого эффекта: как в случае с сетью ВКонтакте: мы могли бы отсеивать тех, кто указывал свою сеть ВКонтакте как основную, и тем самым теряли бы специалистов.
Эта же логика может быть заложена в указанном посте "Психологические различия дам в официальном и гражданском браке" - я показываю, что среди дам в официальном и гражданском браках есть значимые различия по качеству интроверсия - экстраверсия. Но не факт, что это действительно так. Интроверсия - экстраверсия коррелирует с возрастом, а в гражданском браке могут быть просто более молодые женщины, поэтому различие в интроверсия - экстраверсии может быть вызвано не семейным статусом, а возрастом.
Вот это и надо проверять.

Далее код

Данные берем из пакета HRanalitycs 

library(HRanalytics)
data("survey")

Сабсетинг я не буду показывать, если вам, как спецу интересен этот пост, значит вы уже умеете делать это. Я беру полные строки по полу (обратите внимание, там пол указан в разных регистрах), выбираю семейный статус женат и гражданский брак, и полные строки по году рождения. Всего 541 строка. Репрезентативно? Проверим наличие корреляции между возрастом и инстроверсией

cor.test(q1$Год.Рождения, q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.)
 Pearson's product-moment correlation

data:  q1$Год.Рождения and q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.
t = 4.0552, df = 539, p-value = 5.749e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.08905551 0.25269745
sample estimates:
      cor 
0.1720633 

Действительно, слабая (< 0, 3) связь есть. Чем моложе женщина, тем экстравертированней.

Может возникнуть сомнение в нормальности распределения: skew kurtosis для возраста -0.29    -0.12, для интроверсии 0.01    -0.44.
Я знаю три варианта проверки эффекта модерации:

1. в уравнении регрессии предиктор теряем значимость в присутствии модератора
2. формула sobel()
3. SEM - structural equation modeling

Последним я никогда не занимался, а вот первые два метода мы применим. Проблема в том, что у нас аутпут - категориальная переменная, и я честно признаюсь, не знаю, если специальные уравнения для проверки логистической регрессии, Но тем не менее

summary(glm(q1$Семейное.положение.на.момент.трудоустройства ~ q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.+ q1$year, family= 'binomial'))
q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.  -2.0880     0.5396  -3.870 0.000109 ***
q1$year                              -2.8793     0.6573  -4.380 1.18e-05 ***

Показываю уравнение лог регрессии и значимость предикторов. Сами переменные шкалированы в границах 0:1. Значимы обе переменные. Тем не менее, применим формулу sobel(), где зависимая переменная будет as.integer(q1$Семейное.положение.на.момент.трудоустройства)

library(multilevel)
sobel(pred=q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.,med=q1$year,out=q1$Семейное.положение.на.момент.трудоустройства)

У нас q1$Семейное.положение.на.момент.трудоустройства) - это output или зависимая переменная
интроверсия - проверяемая переменная
year - год рождения - медиатор.
Результаты

$`Mod1: Y~X`
             Estimate Std. Error   t value      Pr(>|t|)
(Intercept)  1.967276 0.05170370 38.049047 8.300863e-155
pred        -0.432824 0.09323694 -4.642194  4.334812e-06

$`Mod2: Y~X+M`
              Estimate Std. Error   t value      Pr(>|t|)
(Intercept)  2.2387412 0.07973573 28.077014 1.725524e-107
pred        -0.3620601 0.09306240 -3.890509  1.125716e-04
med         -0.4763393 0.10778717 -4.419258  1.198319e-05

$`Mod3: M~X`
             Estimate Std. Error   t value      Pr(>|t|)
(Intercept) 0.5698981 0.02031517 28.052834 1.948845e-107
pred        0.1485578 0.03663422  4.055165  5.749226e-05

$Indirect.Effect
[1] -0.07076391

Первое, на что обращаем внимание
$Indirect.Effect - сравниваем его со значением 0, 05 - чуть чуть не дотягивает до того, чтобы мы признали наличие эффекта модерации. А счастье было так близко.
Обратите внимание, в уравнении $`Mod2: Y~X+M` значимы также обе переменные.
Поэтому мы делаем, что интроверсия связана с семейным статусом женщины: более интравертированные женщины склонны к официальному браку, экстраверты - к гражданскому.

Здесь

• по оси X - год рожления
• по оси Y - интроверсия - экстраверсия (ближе к десяти - экстраверты)
• цвет точек - семейное положение, красные гражданский брак, зеленые - официальный. 

Наклон линии скорее говорит об эффекте модерации, но это уже другая история, может быть потом напишу еще один пост, и эффекта модерации здесь нет:)