Interaction Effects

Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences by Jacob Cohen, Patricia Cohen, Stephen G. West, Leona S. Aiken

Author

Sungkyun Cho

Published

June 7, 2024

기존 선형모형에서는 각 예측변수들이 additive하게 종속변수에 영향을 미친다고 가정함.

\(\hat{Y} = b_0 + b_1X_1 + b_2X_2 + \cdots + b_kX_k\)

만약, 두 개 이상이 예측변수들이 일종의 시너지 효과와 같이 additive한 효과를 넘어서서 작용한다면,
그 효과를 잡아낼 수 있도록 모형을 변형해야 함.

Continuous vs. Continuous

예제 1

나이가 듦(age)에 따른 지구력(endurance)의 감소가 운동을 한 기간(exercise)에 따라 변화하는가? (p.275)

endurance: the number of minutes of sustained jogging on a treadmill
exercise: the number of years of vigorous physical exercise

  • 연구자의 관심변수에 따라 다르게 표현될 수 있음.
  • 나이가 지구력에 미치는 부정적 영향을 운동이 완화시키는지 관심; 보효요인 (protective factor) <-> 위험요인 (risk factor)
  • 이 때, 운동 기간을 moderator (조절변수)라고 말하고, 그 moderating effect (조절효과)를 가지는지 검증.
  • 통계적으로는 나이과 운동기간이 서로 상호작용(interact)하여 지구력에 영향을 미치는 것으로 나타남.

Data: c07e01dt

acad2 <- read_csv('cohen/data/c07e01dt.csv')
acad2
lowerCor(acad2)  # library(psych)
# A tibble: 245 × 3
     age exercise endurance
   <dbl>    <dbl>     <dbl>
 1    60       10        18
 2    40        9        36
 3    29        2        51
 4    47       10        18
 5    48        9        23
 6    42        6        30
 7    55        8         8
 8    43       19        40
 9    39        9        28
10    51       14        15
# ℹ 235 more rows
          age   exrcs endrn
age        1.00            
exercise   0.28  1.00      
endurance -0.13  0.34  1.00
acad2 %>% 
  ggplot(aes(x = age, y = endurance)) +
  geom_point() +
  geom_smooth(se = FALSE, span = 1)
acad2 %>% 
  ggplot(aes(x = exercise, y = endurance)) +
  geom_point() +
  geom_smooth(se = FALSE, span = 1)

나이와 지구력의 관계가 운동기간에 따라 다른가를 보기 위해, 운동기간을 3구간으로 나눔.

# 편의상 운동기간을 3구간으로 나눔
acad2 %>% 
  mutate(exercise_cat = cut_number(exercise, 3)) %>% 
  ggplot(aes(x = age, y = endurance, color = exercise_cat)) +
  geom_point() +
  geom_smooth(se = FALSE, span = 1) +
  facet_wrap(~exercise_cat)

# 편의상 나이를 3구간으로 나눔
acad2 %>% 
  mutate(age_cat = cut_number(age, 3)) %>% 
  ggplot(aes(x = exercise, y = endurance, color = age_cat)) +
  geom_point() +
  geom_smooth(se = FALSE, span = 1) +
  facet_wrap(~age_cat)

두 변수, 나이와 운동기간으로 지구력을 예측하는 모형을 세우면,

mod_s1 <- lm(endurance ~ age, data = acad2)
mod_s2 <- lm(endurance ~ exercise, data = acad2)
mod <- lm(endurance ~ age + exercise, data = acad2)
export_summs(mod_s1, mod_s2, mod)
Model 1Model 2Model 3
(Intercept)33.16 ***18.39 ***29.40 ***
(3.42)   (1.60)   (3.21)   
age-0.13 *         -0.26 ***
(0.07)          (0.07)   
exercise       0.76 ***0.92 ***
       (0.14)   (0.14)   
N245       245       245       
R20.02    0.11    0.17    
*** p < 0.001; ** p < 0.01; * p < 0.05.

Interaction term을 추가해서 모형을 세우면

\(\displaystyle \widehat{endurance} = b_1\cdot age + b_2\cdot exercise + b_3\cdot age \cdot exercise + b_0\)
      \(\displaystyle = (b_1 + b_3 \cdot exercise)\cdot age + b_2\cdot exercise + b_0\)
이제 age의 기울기가 exercise의 값에 따라 변할 수 있음.

mod_interact <- lm(endurance ~ age * exercise, data = acad2)  
# 동일: endurance ~ age + exercise + age:exercise
S(mod_interact) # library(car)
Call: lm(formula = endurance ~ age * exercise, data = acad2)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  53.17896    7.52661   7.065 1.71e-11 ***
age          -0.76596    0.15980  -4.793 2.87e-06 ***
exercise     -1.35095    0.66626  -2.028 0.043694 *  
age:exercise  0.04724    0.01359   3.476 0.000604 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 9.7 on 241 degrees of freedom
Multiple R-squared: 0.2061
F-statistic: 20.86 on 3 and 241 DF,  p-value: 4.764e-12 
    AIC     BIC 
1814.57 1832.07 

\(\displaystyle \widehat{endurance} = -0.766\cdot age -1.350\cdot exercise + 0.047\cdot age \cdot exercise + 53.18\)
      \(\displaystyle = (-0.766 + 0.047 \cdot exercise)\cdot age -1.350\cdot exercise + 53.18\)

  • age의 기울기는 exercise의 값에 따라 변함
  • 운동기간이 0년인 경우, 지구력에 미치는 나이의 효과: \((-0.766+0.047*0)\cdot age = -0.766\cdot age\)
  • 운동기간이 10년인 경우, 지구력에 미치는 나이의 효과: \((-0.766+0.047*10)\cdot age = -0.30\cdot age\)

각 회귀계수의 의미는 다른 변수의 값이 0일 때의 기울기/효과임

  • age: -0.766은 운동기간이 0년일 때의 나이의 효과
  • exercise: -1.350은 나이가 0살일 때의 운동기간의 효과
  • age:exercise: 0.047는 두 변수의 joint effect
  • 절편 53.18: 운동기간이 0년인 0세의 지구력

회귀계수를 용이하게 해석하기 위해 변수를 centering하거나 standardizing하는 것이 좋음.

  • 0의 의미가 있는 특별한 경우가 아니면 centering을 기본적으로 함.
  • 예를 들어, 언어발달 ~ 나이 * 형제자매 수 + 부모의 교육수준
    • 형제자매의 수 0은 의미가 있음.
  • 위의 경우 ageexercise의 단위가 의미가 있으므로, 표준화보다는 centering
  • p-value들은 모두 바뀐 의미의 회귀계수에 대한 영가설 검정
Important

기존의 additive한 모형에서의 해석을 그대로 적용하면 안됨.

\(\displaystyle \widehat{endurance} = b_0 + b_1\cdot age + b_2 \cdot exercise + b_3 \cdot age \cdot exercise\)

  • age의 효과는 exercise가 고정되었을 때의 효과가 아님! (hold it constant, control for it)
  • 다른 변수를 특정 값에 고정했을 때의 조건부 효과라고 해석할 수 있음.
  • 여기서 \(b_1\)exercise가 0일 때의 age의 효과라고 말할 수 있음.
# jtools의 center()함수를 이용하거나 데이터셋에서 미리 변환
mod_interact_c <- lm(endurance ~ center(age) * center(exercise), data = acad2)
S(mod_interact_c)
Call: lm(formula = endurance ~ center(age) * center(exercise), data = acad2)

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  25.88872    0.64662  40.037  < 2e-16 ***
center(age)                  -0.26169    0.06406  -4.085 6.01e-05 ***
center(exercise)              0.97272    0.13653   7.124 1.20e-11 ***
center(age):center(exercise)  0.04724    0.01359   3.476 0.000604 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 9.7 on 241 degrees of freedom
Multiple R-squared: 0.2061
F-statistic: 20.86 on 3 and 241 DF,  p-value: 4.764e-12 
    AIC     BIC 
1814.57 1832.07 
  • 나이에 대한 회귀계수 -0.26: 운동기간이 평균(10.67년)일 때의 나이의 효과
  • 운동기간에 대한 회귀계수 0.97: 나이가 평균(49.18세)일 때의 운동기간의 효과
psych::describe(acad2)
          vars   n  mean    sd median trimmed   mad min max range skew kurtosis
age          1 245 49.18 10.11     48   49.11 10.38  20  82    62 0.15    -0.08
exercise     2 245 10.67  4.78     11   10.56  4.45   0  26    26 0.27     0.23
endurance    3 245 26.53 10.82     27   26.39 10.38   0  55    55 0.11    -0.30
            se
age       0.65
exercise  0.31
endurance 0.69
library(effects)
# plot(predictorEffect("age", mod_interact, xlevels = 3)) 또는
plot(predictorEffects(mod_interact, ~age, xlevels = 3))  # centering하지 않은 모형

확인하고자 하는 값들을 대입할 수 있음; 보통 평균과에서 +1, -1 표준편차 만큼의 값을 대입함.

m <- mean(acad2$exercise, na.rm = TRUE) |> round(2)
sd <- sd(acad2$exercise, na.rm = TRUE) |> round(2)

plot(predictorEffects(mod_interact, ~age, xlevels = list(exercise = c(m-sd, m, m+sd))))

평균보다 1 표준편차 만큼 오래 운동한 경우: m + 1sd = 10.67 + 4.78 = 15.45년

  • 나이에 따른 지구력 감소가 거의 없는가?
mod_interact_t <- lm(endurance ~ age * I(exercise - 15.45), data = acad2)
S(mod_interact_t)
Call: lm(formula = endurance ~ age * I(exercise - 15.45), data = acad2)

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)             32.30671    4.72744   6.834 6.71e-11 ***
age                     -0.03602    0.09027  -0.399 0.690184    
I(exercise - 15.45)     -1.35095    0.66626  -2.028 0.043694 *  
age:I(exercise - 15.45)  0.04724    0.01359   3.476 0.000604 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 9.7 on 241 degrees of freedom
Multiple R-squared: 0.2061
F-statistic: 20.86 on 3 and 241 DF,  p-value: 4.764e-12 
    AIC     BIC 
1814.57 1832.07 

가설 검정(p-value)보다는 신뢰구간(confidence interval)을 보는 것을 항상 권장함.

  • R의 confint()나 car 패키지의 Confint() 함수를 이용
  • jtools에서 summ(mod, confint = TRUE)
Confint(mod_interact_t) |> print(digits = 2)  # car 패키지
                        Estimate 2.5 % 97.5 %
(Intercept)               32.307 22.99 41.619
age                       -0.036 -0.21  0.142
I(exercise - 15.45)       -1.351 -2.66 -0.039
age:I(exercise - 15.45)    0.047  0.02  0.074

모든 moderator의 값에 대한 통계적 유의도를 계산

# Johnson-Neyman intervals 
library(interactions)
johnson_neyman(mod_interact, pred = "age", modx = "exercise")
JOHNSON-NEYMAN INTERVAL 

When exercise is OUTSIDE the interval [13.21, 24.37], the slope of age is p
< .05.

Note: The range of observed values of exercise is [0.00, 26.00]

회귀모형에 대한 진단

avPlots(mod_interact_c, id=list(n=3))

mod_interact_update <- update(mod_interact_c, subset = -c(3, 34)) # 3, 34번째 관측치 제거
avPlots(mod_interact_update)

S(mod_interact_update)
Call: lm(formula = endurance ~ center(age) * center(exercise), data = acad2,
         subset = -c(3, 34))

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  25.88331    0.63738  40.609  < 2e-16 ***
center(age)                  -0.25750    0.06371  -4.041 7.17e-05 ***
center(exercise)              0.95987    0.13561   7.078 1.61e-11 ***
center(age):center(exercise)  0.03482    0.01402   2.483   0.0137 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 9.561 on 239 degrees of freedom
Multiple R-squared: 0.195
F-statistic:  19.3 on 3 and 239 DF,  p-value: 3.063e-11 
    AIC     BIC 
1792.81 1810.27 

상호작용 효과의 크기

  • 상호작용 항이 추가된 모형과 추가되지 않은 모형을 비교
  • \(\Delta R^2 = R^2_{full} - R^2_{reduced} = 0.21 - 0.17 = 0.04\)
  • 상호작용이 추가된 모형이 4%의 추가적인 설명력을 가짐.
  • 이 크기에 대한 가설 검정은 상호작용항의 회귀계수에 대한 검정과 동일; 함수 anova()이용
mod <- lm(endurance ~ age + exercise, data = acad2)
mod_interact <- lm(endurance ~ age * exercise, data = acad2)
export_summs(mod, mod_interact, error_format = "({p.value})", number_format = "%.4f")
Model 1Model 2
(Intercept)29.3952 ***53.1790 ***
(0.0000)   (0.0000)   
age-0.2571 ***-0.7660 ***
(0.0001)   (0.0000)   
exercise0.9163 ***-1.3510 *  
(0.0000)   (0.0437)   
age:exercise         0.0472 ***
         (0.0006)   
N245         245         
R20.1663    0.2061    
*** p < 0.001; ** p < 0.01; * p < 0.05.
anova(mod, mod_interact)
Analysis of Variance Table

Model 1: endurance ~ age + exercise
Model 2: endurance ~ age * exercise
  Res.Df   RSS Df Sum of Sq     F    Pr(>F)    
1    242 23810                                 
2    241 22674  1    1136.5 12.08 0.0006042 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model 1

  • 디폴트로 moderator의 16th, 50th, 84th 세 값에 대한 조건부 효과를 계산
  • moments=1: moderator의 mean, mean+sd, mean-sd 세 값에 대한 조건부 효과를 계산
  • wmodval=c(a, b, c)는 moderator 값이 a, b, c일 때의 조건부 효과를 계산
  • jn=1: Johnson-Neyman interval 계산
  • center=1: centering
process(data=acad2, y="endurance", x="age", w="exercise", model=1, center=1, jn=1)

# ********************* PROCESS for R Version 4.3.1 ********************* 
 
#            Written by Andrew F. Hayes, Ph.D.  www.afhayes.com              
#    Documentation available in Hayes (2022). www.guilford.com/p/hayes3   
 
# *********************************************************************** 
                 
# Model : 1        
#     Y : endurance
#     X : age      
#     W : exercise 

# Sample size: 245


# *********************************************************************** 
# Outcome Variable: endurance

# Model Summary: 
#           R      R-sq       MSE         F       df1       df2         p
#      0.4540    0.2061   94.0821   20.8585    3.0000  241.0000    0.0000

# Model: 
#              coeff        se         t         p      LLCI      ULCI
# constant   25.8887    0.6466   40.0371    0.0000   24.6150   27.1625
# age        -0.2617    0.0641   -4.0848    0.0001   -0.3879   -0.1355
# exercise    0.9727    0.1365    7.1244    0.0000    0.7038    1.2417
# Int_1       0.0472    0.0136    3.4757    0.0006    0.0205    0.0740

# Product terms key:
# Int_1  :  age  x  exercise      

# Test(s) of highest order unconditional interaction(s):
#       R2-chng         F       df1       df2         p
# X*W    0.0398   12.0804    1.0000  241.0000    0.0006
# ----------
# Focal predictor: age (X)
#       Moderator: exercise (W)

# Conditional effects of the focal predictor at values of the moderator(s):
#    exercise    effect        se         t         p      LLCI      ULCI
#     -4.6735   -0.4825    0.0911   -5.2935    0.0000   -0.6620   -0.3029
#      0.3265   -0.2463    0.0641   -3.8403    0.0002   -0.3726   -0.1199
#      4.3265   -0.0573    0.0861   -0.6656    0.5063   -0.2268    0.1123

# Moderator value(s) defining Johnson-Neyman significance region(s):
#       Value   % below   % above
#      2.5328   73.8776   26.1224
#     13.6954   99.1837    0.8163

# Conditional effect of focal predictor at values of the moderator:
#    exercise    effect        se         t         p      LLCI      ULCI
#    -10.6735   -0.7660    0.1598   -4.7931    0.0000   -1.0807   -0.4512
#     -9.3050   -0.7013    0.1430   -4.9057    0.0000   -0.9829   -0.4197
#     -7.9366   -0.6367    0.1266   -5.0288    0.0000   -0.8860   -0.3873
#     -6.5682   -0.5720    0.1110   -5.1552    0.0000   -0.7906   -0.3534
#     -5.1998   -0.5074    0.0964   -5.2647    0.0000   -0.6972   -0.3175
#     -3.8314   -0.4427    0.0834   -5.3087    0.0000   -0.6070   -0.2784
#     -2.4629   -0.3781    0.0729   -5.1863    0.0000   -0.5216   -0.2345
#     -1.0945   -0.3134    0.0661   -4.7437    0.0000   -0.4435   -0.1833
#      0.2739   -0.2487    0.0641   -3.8809    0.0001   -0.3750   -0.1225
#      1.6423   -0.1841    0.0674   -2.7312    0.0068   -0.3169   -0.0513
#      2.5328   -0.1420    0.0721   -1.9699    0.0500   -0.2841    0.0000
#      3.0107   -0.1194    0.0753   -1.5862    0.1140   -0.2678    0.0289
#      4.3792   -0.0548    0.0865   -0.6332    0.5272   -0.2253    0.1157
#      5.7476    0.0099    0.1000    0.0985    0.9216   -0.1871    0.2069
#      7.1160    0.0745    0.1149    0.6484    0.5174   -0.1519    0.3009
#      8.4844    0.1392    0.1308    1.0641    0.2883   -0.1184    0.3967
#      9.8528    0.2038    0.1473    1.3839    0.1677   -0.0863    0.4939
#     11.2213    0.2685    0.1642    1.6348    0.1034   -0.0550    0.5919
#     12.5897    0.3331    0.1815    1.8354    0.0677   -0.0244    0.6906
#     13.6954    0.3853    0.1956    1.9699    0.0500    0.0000    0.7707
#     13.9581    0.3978    0.1990    1.9988    0.0468    0.0058    0.7898
#     15.3265    0.4624    0.2167    2.1339    0.0339    0.0356    0.8893

# ******************** ANALYSIS NOTES AND ERRORS ************************ 

# Level of confidence for all confidence intervals in output: 95

# W values in conditional tables are the 16th, 50th, and 84th percentiles.
 
# NOTE: The following variables were mean centered prior to analysis: 
#          exercise age
Interaction의 패턴
  1. Synergistic or enhancing interaction
  • 상호작용 효과가 원래 효과들과 같은 방향으로 작용하는 경우
  • 삶의 만족도(Y)가 직업 스트레스(X)와 부정적인 관계에 있고, 부부관계의 문제(Z)와도 부정적인 관계에 있는 경우
  • 이 둘의 상호작용이 부정적이라면, 직업 스트레스와 부부관계의 문제가 동시에 증가하면 각각의 sum이 예측하는 것보다 더 낮은 삶의 만족도가 예측됨.
  1. Buffering interaction
  • 두 변수가 반대 방향으로 Y에 작용하고 있을 때, 한 변수가 다른 변수의 효과를 감소시키는 경우
  • 즉, 한 변수의 impact가 다른 변수의 impact를 줄여주는 경우
  • 건강보건에 대한 연구에서, 한 변수가 질병의 위험요인이고 다른 변수가 질병의 위험을 줄여주는 보호요인인 경우
  • 위의 예에서처럼, 나이(X)는 지구력 감소의 위험요인이고, 운동기간(Z)은 지구력 보호요인인 경우
  1. Interference or antagonistic interactionin
  • 두 변수가 같은 방향으로 Y에 작용하고 있을 때, 상호작용은 반대 방향으로 작용하는 경우
  • 대학생의 학업성취도(Y)에 대하여, 학업동기(X)와 학업능력(Z)이 모두 학업성취도(Y)에 긍정적인 영향을 미치나 이 두 변수는 서로 보완적인 효과를 가지고 있음.
  • 즉, 성취도에 대한 학업능력의 중요성은 높은 학업동기에 의해 낮아질 수 있음.
  • 반대로, 학업동기에 대한 중요성은 높은 학업능력에 의해 낮아질 수 있음.

예제 2

Source: p.551, Statistical Methods for Psychology (8e) by Dave C. Howell

오리엔테이션에 참석한 대학 신입생을 대상으로, 스트레스 받은 일(hassles)이 많을 수록 여러 증상들(symptoms)을 더 경험하는데, 주위의 지지(support)가 많을 수록 그 증상들이 감소한다는 것을 알아보고자 설문 조사. (Wagner, Compas, and Howell (1988))

  • hassles: 소소한 일상적인 스트레스를 경험한 횟수
  • support: 본인이 인지하는 주위의 지지 정도
  • symptoms: 증상에 대한 checklist로 나타난 증상의 개수

Data: hassles.csv

hassles <- read_csv("data/hassles.csv")
hassles
# A tibble: 56 × 3
   hassles support symptoms
     <dbl>   <dbl>    <dbl>
 1     176      10       73
 2     379      50       88
 3     126      45      118
 4     193      40       79
 5     229      40      127
 6     153      39       73
 7     214      38       93
 8     164      37       99
 9     143      36       81
10      27      36       64
# ℹ 46 more rows
car::scatterplotMatrix(hassles)

hassles |> 
  lowerCor()  # library(psych)
         hssls spprt sympt
hassles   1.00            
support  -0.17  1.00      
symptoms  0.58 -0.13  1.00
hassles <- hassles |> 
  mutate(
    hassles_c = center(hassles),  # library(jtools)
    support_c = center(support)
  )

\(\hat{symptom} = b_0 + b_1\cdot hassles_c + b_2\cdot support_c + b_3\cdot hassles_c \cdot support_c\)

mod <- lm(symptoms ~ hassles_c * support_c, data = hassles)
S(mod)  # library(car)
Call: lm(formula = symptoms ~ hassles_c * support_c, data = hassles)

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         89.584940   2.291503  39.094  < 2e-16 ***
hassles_c            0.085942   0.019213   4.473 4.22e-05 ***
support_c            0.146358   0.305244   0.479   0.6336    
hassles_c:support_c -0.005065   0.002363  -2.144   0.0368 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 16.89 on 52 degrees of freedom
Multiple R-squared: 0.3885
F-statistic: 11.01 on 3 and 52 DF,  p-value: 1.046e-05 
   AIC    BIC 
481.39 491.51 
hassles |> ggplot(aes(x = support_c)) + geom_freqpoly()

library(effects)
# 임의의 3 levels 선택
predictorEffects(mod, xlevels = 3) |>  # 임의의 3 levels 선택
  plot(lines = list(multiline = TRUE)) # 3개의 선을 함께 그림

# 특정 3 levels 선택: -10(평균-10), 0(평균), 10(평균+10)
predictorEffects(mod, ~ hassles_c, 
                 xlevels = list(support_c = c(-10, 0, 10))) |>  # 특정 3 levels 선택
  plot(lines = list(multiline = TRUE))  # 3개의 선을 함께 그림

support_c가 10일 때, hassels의 기울기는 통계적으로 유의한가를 테스트

mod_probe <- lm(symptoms ~ hassles_c * I(support_c - 10), data = hassles)
S(mod_probe)  # library(car)
Call: lm(formula = symptoms ~ hassles_c * I(support_c - 10), data = hassles)

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 91.048515   3.697199  24.626   <2e-16 ***
hassles_c                    0.035293   0.034034   1.037   0.3045    
I(support_c - 10)            0.146358   0.305244   0.479   0.6336    
hassles_c:I(support_c - 10) -0.005065   0.002363  -2.144   0.0368 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 16.89 on 52 degrees of freedom
Multiple R-squared: 0.3885
F-statistic: 11.01 on 3 and 52 DF,  p-value: 1.046e-05 
   AIC    BIC 
481.39 491.51 

모든 moderator의 값에 대한 통계적 유의도를 계산

# Johnson-Neyman intervals 
library(interactions)
johnson_neyman(mod, pred = "hassles_c", modx = "support_c")
JOHNSON-NEYMAN INTERVAL 

When support_c is OUTSIDE the interval [6.25, 297.50], the slope of
hassles_c is p < .05.

Note: The range of observed values of support_c is [-18.96, 21.04]

회귀모형에 대한 진단

avPlots(mod)

influenceIndexPlot(mod)

mod_update <- update(mod, subset = -52) # 52번째 관측치 제거
avPlots(mod_update)

S(mod_update)
Call: lm(formula = symptoms ~ hassles_c * support_c, data = hassles, subset =
         -52)

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         89.140023   2.156778  41.330  < 2e-16 ***
hassles_c            0.070908   0.018801   3.772 0.000423 ***
support_c            0.125998   0.286624   0.440 0.662089    
hassles_c:support_c -0.001655   0.002524  -0.656 0.515015    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 15.86 on 51 degrees of freedom
Multiple R-squared: 0.2291
F-statistic: 5.053 on 3 and 51 DF,  p-value: 0.00385 
   AIC    BIC 
465.93 475.97 
Linear vs. Nonlinear Interaction

위에서 다룬 예제들은 모두 선형적인 상호작용을 다룸; 가장 단순한 형태의 상호작용
다음과 같은 비선형적인 상호작용도 가능함. (p. 292)

\(X\)에 대한 2차 다항함수로 fit을 한다면,
이 때, 상호작용이 없다면
\(\hat{Y} = b_0 + b_1X + b_2X^2 + b_3Z\)

상호작용을 추가하면
\(\hat{Y} = b_0 + b_1X + b_2X^2 + b_3XZ + b_4X^2Z\)

Formula: Y ~ X + I(X^2) + Z + X:Z + I(X^2):Z

Continuous vs. Categorical

  • 성별 이나 실험조건/통제조건 같이 두 개의 범주를 가지는 변수와의 상호작용
  • 직업군과 같이 세 개 이상의 범주를 가지는 변수와의 상호작용

Binary 변수와의 상호작용

예제 1: Data from the 1985 Current Population Survey (CPS85)
나이(age)에 따른 임금(wage)의 증가량이 남녀(sex)에 따라 다른가?

cps <- mosaicData::CPS85 |> as_tibble()
cps |> head(3)
# A tibble: 3 × 11
   wage  educ race  sex   hispanic south married exper union   age sector
  <dbl> <int> <fct> <fct> <fct>    <fct> <fct>   <int> <fct> <int> <fct> 
1   9      10 W     M     NH       NS    Married    27 Not      43 const 
2   5.5    12 W     M     NH       NS    Married    20 Not      38 sales 
3   3.8    12 W     F     NH       NS    Single      4 Not      22 sales 

편의상 management 섹터에서만 보면,

cps <- cps |> 
  filter(sector == "manag") |> 
  filter(wage < 30 & age < 60) |> 
  slice(-51)
cps |> ggplot(aes(x = age, y = wage, color = sex)) + geom_point() + geom_smooth(se=F)

\(\hat{wage} = b_0 + b_1\cdot age + b_2\cdot sexM + b_3\cdot age \cdot sexM\)

mod <- lm(wage ~ age * sex, data = cps)
S(mod)
Call: lm(formula = wage ~ age * sex, data = cps)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   4.1600     3.9354   1.057   0.2960  
age           0.1526     0.1044   1.463   0.1504  
sexM         -6.7728     5.4323  -1.247   0.2188  
age:sexM      0.2852     0.1409   2.025   0.0487 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 4.795 on 46 degrees of freedom
Multiple R-squared: 0.4266
F-statistic: 11.41 on 3 and 46 DF,  p-value: 1.028e-05 
   AIC    BIC 
304.47 314.03 
predictorEffects(mod) |> plot()

상호작용항의 추가 설명력: \(\Delta R^2\) 계산

mod <- lm(wage ~ age * sex, data = cps)
mod_reduced <- lm(wage ~ age + sex, data = cps)

r2_full <- summary(mod)$r.squared
r2_reduced <- summary(mod_reduced)$r.squared

sprintf("R2 diff: %.3f", r2_full - r2_reduced)
[1] "R2 diff: 0.051"
anova(mod_reduced, mod)
Analysis of Variance Table

Model 1: wage ~ age + sex
Model 2: wage ~ age * sex
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     47 1151.7                              
2     46 1057.5  1     94.23 4.0991 0.04874 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

예제 2: 기후 변화 회의론이 원조 보류의 정당성에 미치는 영향

p. 285, Introduction to Mediation, Moderation, and Conditional Process Analysis (3e) by Andrew F. Hayes

저는 7장에서 소개한 조절 분석 방법을 재난의 원인을 기후 변화로 돌리는 것과 원인을 특정하지 않는 것이 기후 변화의 현실에 대한 사람들의 믿음에 따라 피해자에 대한 지원을 보류하는 사람들의 정당성에 차별적으로 영향을 미치는지 조사한 연구 데이터를 사용하여 설명했습니다. 이 예에서 초점 선행 변수는 가뭄의 원인에 대한 프레임을 코딩하는 이분법적 변수였고, 조절 변수는 기후 변화의 현실에 대한 회의론의 연속선상에 각 사람을 위치시키는 측정된 개인적 차이였습니다. 분석 결과, 기후 변화에 대해 회의적인 사람들은 원인을 특정하지 않은 경우보다 기후 변화로 인한 재난일 때 원조 보류에 대한 더 강력한 정당성을 보고했습니다. 기후 변화에 대한 회의적 시각이 상대적으로 낮은 사람들 사이에서는 이러한 효과가 관찰되지 않았습니다. 하지만 재난의 원인을 어떻게 규정하느냐가 아니라 기후변화 회의론이 재난 피해자에 대한 원조 의지에 미치는 영향에 실질적으로 초점을 맞춘다면 어떨까요? 이 질문은 간단한 회귀 분석으로 쉽게 답할 수 있습니다. 기후 변화 회의론(X)으로부터 원조 보류의 정당성(Y)을 추정하는 가장 적합한 OLS 회귀 모델은 Yˆ = 2.186+0.201X입니다. 따라서 기후 변화에 대한 회의론이 1단위 차이가 나는 두 사람은 기근 피해자에 대한 원조 보류의 정당성 강도에서 0.201단위 차이가 나는 것으로 추정됩니다. 이 관계는 통계적으로 유의미합니다. 그러나 이 분석은 이 연구 참여자의 절반은 가뭄과 그로 인한 기근이 기후 변화로 인한 것이라고 들었지만 다른 사람들은 그 원인에 대한 정보가 전혀 제공되지 않았다는 점을 완전히 무시하고 있습니다. 이 연구의 저자에 따르면, 기후 변화로 인한 것으로 표시된 사건은 기후 변화 회의론자들이 기후 변화의 현실에 대해 덜 회의적인 사람들에 비해 피해자를 돕는 것의 가치에 대해 특히 의심하게 만들 수 있다고 합니다. 즉, 기후 변화에 대한 귀인이 동기가 부여된 방어적 태도를 유발할 수 있으며, 기후 변화의 현실에 대한 태도는 그러한 귀인이 없을 때와 비교하여 도움에 대한 태도를 더 잘 예측할 수 있습니다. 다시 말해, 이러한 귀인은 사람들이 아무런 원인이 제공되지 않았을 때보다 자신의 태도와 더 일치하는 방식으로 피해자의 요구에 응답하도록 유도할 수 있습니다. 이러한 추론에 따르면 기후변화가 가뭄의 원인이라고 들었을 때와 원인에 대한 정보가 제공되지 않았을 때 기후변화 회의론과 원조 보류의 정당성 사이의 관계가 다를 것으로 예상할 수 있습니다. 기후변화 회의론 X와 재난 프레임 W(기후변화 조건은 1, 자연적 원인 조건은 0으로 설정)라고 부르면 7장에서와 마찬가지로 간단한 조정 모델을 추정할 수 있는데, 이 모델에서 X가 Y에 미치는 영향은 다음과 같다. 이러한 과정은 그림 8.1, 패널 A에 개념적으로 도식화되어 있으며 그림 8.1, 패널 B의 통계 다이어그램에서와 같이 X, W 및 XW를 선행 변수로 하는 통계 모형으로 변환됩니다.

Translated with DeepL.com (free version)

disaster <- read_csv("data/hayes2022data/disaster/disaster.csv")
disaster <- disaster |> 
  mutate(frame_f = factor(frame, levels = c(0, 1), labels = c("control", "framed")))
disaster |> head(3)
# A tibble: 3 × 6
     id frame donate justify skeptic frame_f
  <dbl> <dbl>  <dbl>   <dbl>   <dbl> <fct>  
1     1     1    5.6    2.95     1.8 framed 
2     2     1    4.2    2.85     5.2 framed 
3     3     1    4.2    3        3.2 framed 
Show the code
disaster |> 
  ggplot(aes(x = skeptic, y = justify, color = frame_f)) +
  geom_point() +
  geom_smooth(se = F, method = "lm", color="grey50") +
  geom_smooth(se = F) +
  facet_wrap(~frame_f)

\(\hat{justify} = b_0 + b_1\cdot skeptic + b_2\cdot frame + b_3\cdot skeptic \cdot frame\)

mod <- lm(justify ~ center(skeptic) * frame_f, data = disaster)
S(mod)
Call: lm(formula = justify ~ center(skeptic) * frame_f, data = disaster)

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    2.80650    0.07753  36.201  < 2e-16 ***
center(skeptic)                0.10508    0.03813   2.756 0.006375 ** 
frame_fframed                  0.11712    0.11206   1.045 0.297156    
center(skeptic):frame_fframed  0.20118    0.05527   3.640 0.000344 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 0.8129 on 207 degrees of freedom
Multiple R-squared: 0.2463
F-statistic: 22.54 on 3 and 207 DF,  p-value: 1.138e-12 
   AIC    BIC 
517.36 534.12 
m <- mean(disaster$skeptic, na.rm = TRUE) |> round(2)
sd <- sd(disaster$skeptic, na.rm = TRUE) |> round(2)

predictorEffects(mod, xlevels = list(skeptic = c(m-sd, m, m+sd))) |>
  plot(lines = list(multiline = TRUE)) # 3개의 선을 함께 그림

# johnson_neyman() 함수의 경우 카테고리 변수를 factor가 아닌 숫자로 입력해야 함
mod2 <- lm(justify ~ skeptic * frame, data = disaster)  # frame: 0, 1
johnson_neyman(mod2, pred = "frame", modx = "skeptic")
JOHNSON-NEYMAN INTERVAL 

When skeptic is OUTSIDE the interval [1.17, 3.93], the slope of frame is p
< .05.

Note: The range of observed values of skeptic is [1.00, 9.00]

상호작용항의 추가 설명력: \(\Delta R^2\) 계산

mod <- lm(justify ~ skeptic * frame_f, data = disaster)
mod_reduced <- lm(justify ~ skeptic + frame_f, data = disaster)

r2_full <- summary(mod)$r.squared
r2_reduced <- summary(mod_reduced)$r.squared

sprintf("R2 diff: %.3f", r2_full - r2_reduced)
[1] "R2 diff: 0.048"

3개 이상의 범주를 가지는 변수와의 상호작용

Data: c0904dt2.csv

acad3 <- read_csv("data/c0904dt2.csv")
acad3 |> head(5)
# A tibble: 5 × 5
  depart    pub  time salary sex   
  <chr>   <dbl> <dbl>  <dbl> <chr> 
1 anthrop    16     3 56465. female
2 anthrop    25     7 92044. male  
3 anthrop    16     2 48980. female
4 anthrop    24     1 53239. female
5 anthrop    24     8 98948. male  
acad3 <- acad3 |> mutate(pub_c = center(pub), time_c = center(time))
acad3 |> 
  ggplot(aes(x = pub, y = salary, color = depart)) +
  geom_point() + geom_smooth(method=lm) + facet_wrap(~depart)

\(\hat{salary} = b_0 + b_1\cdot pub + b_2\cdot depart_{hist} + b_3\cdot depart_{soc} + b_4\cdot pub \cdot depart_{hist} + b_5\cdot pub \cdot depart_{soc}\)

mod <- lm(salary ~ pub_c * depart, data = acad3)
S(mod)
Call: lm(formula = salary ~ pub_c * depart, data = acad3)

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       56917.9     2207.4  25.785  < 2e-16 ***
pub_c              1372.9      252.5   5.438 2.25e-07 ***
departhist         9796.1     3615.1   2.710  0.00755 ** 
departsoc          9672.8     3235.2   2.990  0.00328 ** 
pub_c:departhist   -961.0      466.2  -2.061  0.04109 *  
pub_c:departsoc   -1115.0      495.4  -2.251  0.02592 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 15670 on 144 degrees of freedom
Multiple R-squared: 0.1892
F-statistic: 6.722 on 5 and 144 DF,  p-value: 1.167e-05 
    AIC     BIC 
3331.44 3352.51 
plot(predictorEffects(mod, ~depart, xlevels = 3))

plot(predictorEffects(mod, ~pub_c))

상호작용항의 추가 설명력: \(\Delta R^2\) 계산

mod <- lm(salary ~ pub_c * depart, data = acad3)
mod_reduced <- lm(salary ~ pub_c + depart, data = acad3)

r2_full <- summary(mod)$r.squared
r2_reduced <- summary(mod_reduced)$r.squared

sprintf("R2 diff: %.3f", r2_full - r2_reduced)
[1] "R2 diff: 0.041"
anova(mod_reduced, mod)
Analysis of Variance Table

Model 1: salary ~ pub_c + depart
Model 2: salary ~ pub_c * depart
  Res.Df        RSS Df  Sum of Sq      F  Pr(>F)  
1    146 3.7161e+10                               
2    144 3.5366e+10  2 1795358985 3.6551 0.02829 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

연차의 효과를 고려한다면

mod2 <- lm(salary ~ (time_c + pub_c) * depart, data = acad3)
S(mod2, brief = T)
Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        59318.9     2222.0  26.696  < 2e-16 ***
time_c              1144.4      387.8   2.951  0.00371 ** 
pub_c               1063.7      258.7   4.112 6.62e-05 ***
departhist          4223.7     3632.6   1.163  0.24691    
departsoc           6999.9     3140.1   2.229  0.02738 *  
time_c:departhist    253.7      598.2   0.424  0.67212    
time_c:departsoc    -151.1      590.0  -0.256  0.79829    
pub_c:departhist    -985.4      462.1  -2.133  0.03470 *  
pub_c:departsoc     -793.6      475.8  -1.668  0.09752 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 14680 on 141 degrees of freedom
Multiple R-squared: 0.3035
F-statistic: 7.679 on 8 and 141 DF,  p-value: 1.655e-08 
    AIC     BIC 
3314.66 3344.77 
mod3 <- lm(salary ~ time_c + pub_c * depart, data = acad3)
S(mod3, brief = T)
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       59377.2     2119.2  28.018  < 2e-16 ***
time_c             1172.2      244.6   4.792 4.09e-06 ***
pub_c              1056.2      244.3   4.324 2.86e-05 ***
departhist         4677.9     3532.7   1.324   0.1876    
departsoc          6892.6     3068.9   2.246   0.0262 *  
pub_c:departhist   -924.0      434.4  -2.127   0.0351 *  
pub_c:departsoc    -783.9      466.6  -1.680   0.0951 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 14600 on 143 degrees of freedom
Multiple R-squared: 0.3014
F-statistic: 10.28 on 6 and 143 DF,  p-value: 1.874e-09 
    AIC     BIC 
3311.10 3335.18 
plot(predictorEffects(mod3, ~depart, xlevels = 3))

Categorical vs. Categorical

ANOVA: factorical design

실험연구의 예

실험연구 분석에 대한 더 논의가 필요함…

Data: C0901DT.csv
2 X 2인 경우: lesion(2 levels) x drug(2 levels)

rats <- read_csv("data/C0901DT.csv")
rats <- rats |> rename_with(str_to_lower)  # 변수명 소문자로 변경
rats |> sample_n(6)
# A tibble: 6 × 5
     ya    yb    yc lesion  drug   
  <dbl> <dbl> <dbl> <chr>   <chr>  
1  5.81 13.8   7.81 SURGERY ACTIVE 
2  9.74  9.74 11.7  SURGERY PLACEBO
3  3.09  3.09  1.09 SHAM    ACTIVE 
4  9.58 13.6  11.6  SHAM    PLACEBO
5  3.86  3.86  1.86 SHAM    ACTIVE 
6  3.92 11.9   5.92 SURGERY ACTIVE 
rats |> group_by(lesion, drug) |> summarise(yc_mean = mean(yc))
# A tibble: 4 × 3
  lesion  drug    yc_mean
  <chr>   <chr>     <dbl>
1 SHAM    ACTIVE     2.00
2 SHAM    PLACEBO    9.99
3 SURGERY ACTIVE     8.02
4 SURGERY PLACEBO   12.0 
mod <- lm(yb ~ lesion * drug, data = rats)
S(mod)
Call: lm(formula = yb ~ lesion * drug, data = rats)

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 3.9965     0.1524   26.23   <2e-16 ***
lesionSURGERY              10.0251     0.2302   43.55   <2e-16 ***
drugPLACEBO                 8.0141     0.2056   38.98   <2e-16 ***
lesionSURGERY:drugPLACEBO -12.0336     0.3256  -36.96   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 0.9758 on 147 degrees of freedom
Multiple R-squared: 0.9401
F-statistic: 768.6 on 3 and 147 DF,  p-value: < 2.2e-16 
   AIC    BIC 
427.07 442.15 
predictorEffects(mod, ~drug) |> plot()

Data: C0902DT.csv
4 X 3인 경우: hospital(4 levels) x treatment(3 levels)

treat <- read_csv("data/C0902DT.csv")
treat <- treat |> rename_with(str_to_lower) # 변수명 소문자로 변경
treat |> head(5)
# A tibble: 5 × 3
  hospital treatmen     y
     <dbl>    <dbl> <dbl>
1        1        1  66.8
2        1        1  42.4
3        1        1  20.7
4        1        1  46.4
5        1        1  54.3
treat <- treat |> 
  mutate(hospital = factor(hospital), treatmen = factor(treatmen))
mod <- lm(y ~ hospital * treatmen, data = treat)
predictorEffects(mod, ~treatmen) |> plot()

관찰연구의 예

성별과 결혼 여부에 따른 임금차이

cps <- mosaicData::CPS85 |> as_tibble()
cps |> head(5)
# A tibble: 5 × 11
   wage  educ race  sex   hispanic south married exper union   age sector  
  <dbl> <int> <fct> <fct> <fct>    <fct> <fct>   <int> <fct> <int> <fct>   
1   9      10 W     M     NH       NS    Married    27 Not      43 const   
2   5.5    12 W     M     NH       NS    Married    20 Not      38 sales   
3   3.8    12 W     F     NH       NS    Single      4 Not      22 sales   
4  10.5    12 W     F     NH       NS    Married    29 Not      47 clerical
5  15      12 W     M     NH       NS    Married    40 Union    58 const   
mod <- lm(wage ~ sex * married, data = cps)
S(mod)
Call: lm(formula = wage ~ sex * married, data = cps)

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)          7.6838     0.3898  19.711  < 2e-16 ***
sexM                 3.1923     0.5319   6.002 3.62e-09 ***
marriedSingle        0.5759     0.6697   0.860  0.39026    
sexM:marriedSingle  -3.0972     0.9073  -3.414  0.00069 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 4.962 on 530 degrees of freedom
Multiple R-squared: 0.07314
F-statistic: 13.94 on 3 and 530 DF,  p-value: 9.234e-09 
    AIC     BIC 
3232.05 3253.45 
predictorEffects(mod) |> plot()

관찰연구의 경우 confounding이 될만한 통제 변수를 충분히 고려해야 하고, 나이가 성별과 결혼 여부와 상관관계가 있으므로 통제 변수로 포함.

cps |> select(age, sex, married) |> lowerCor()
         age   sex*  mrrd*
age       1.00            
sex*     -0.08  1.00      
married* -0.28  0.01  1.00
# age와 wage의 관계는 비선형적: spline regression (piecewise polynomial) fit을 고려할 필요있음
mod2 <- lm(wage ~ sex * married + age + I(age^2), data = cps)  
predictorEffects(mod2, ~ sex + married) |> plot()

S(mod2)
Call: lm(formula = wage ~ sex * married + age + I(age^2), data = cps)

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -5.147951   2.391631  -2.152  0.03181 *  
sexM                3.152425   0.516297   6.106 1.98e-09 ***
marriedSingle       1.321511   0.662608   1.994  0.04662 *  
age                 0.608443   0.121565   5.005 7.62e-07 ***
I(age^2)           -0.006630   0.001483  -4.470 9.60e-06 ***
sexM:marriedSingle -2.617902   0.887868  -2.949  0.00333 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 4.816 on 528 degrees of freedom
Multiple R-squared: 0.1301
F-statistic:  15.8 on 5 and 528 DF,  p-value: 1.649e-14 
    AIC     BIC 
3202.17 3232.13 

결혼 유무에 따른 임금의 차이는 여성의 경우 미혼이 $1.32(p = 0.047) 높음.
남성의 경우는 아래와 같이 factor level을 변경해서 보면, 남성인 경우 미혼이 $1.29(p = 0.041) 낮음.

mod2_m <- lm(wage ~ fct_relevel(sex, "M") * married + age + I(age^2), data = cps)
S(mod2_m)
Call: lm(formula = wage ~ fct_relevel(sex, "M") * married + age + I(age^2),
         data = cps)

Coefficients:
                                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          -1.995525   2.393324  -0.834  0.40478    
fct_relevel(sex, "M")F               -3.152425   0.516297  -6.106 1.98e-09 ***
marriedSingle                        -1.296391   0.632255  -2.050  0.04082 *  
age                                   0.608443   0.121565   5.005 7.62e-07 ***
I(age^2)                             -0.006630   0.001483  -4.470 9.60e-06 ***
fct_relevel(sex, "M")F:marriedSingle  2.617902   0.887868   2.949  0.00333 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 4.816 on 528 degrees of freedom
Multiple R-squared: 0.1301
F-statistic:  15.8 on 5 and 528 DF,  p-value: 1.649e-14 
    AIC     BIC 
3202.17 3232.13 
Important

특히, 카테고리 변수와의 상호작용을 테스트할 때, 유닉크한 카테고리들의 조합에 해당하는 표본의 수가 급격히 줄기 때문에 표본의 수가 충분히 많아야 함.
예를 들어, 위의 경우 4가지 조합(싱글 남자, 싱글 여자, 기혼 남자, 기혼 여자)에 대한 표본의 수가 충분히 많아야 함.

Moderated Mediation

p. 424, Introduction to Mediation, Moderation, and Conditional Process Analysis (3e) by Andrew F. Hayes

11.3 예시: 업무 팀에게 자신의 감정 숨기기 수년에 걸쳐 대중음악은 우리의 직관을 강화해왔고, 친한 친구들이 제공한 조언이나 친한 친구들이 제공하고 토크쇼 심리학자들은 감정을 억누르고 다른 사람의 눈에 띄지 않도록 숨겨서는 좋은 결과를 얻을 수 없다고 강조합니다. 당대의 아티스트들은 “자신을 표현하라”(마돈나)고, “침묵의 강령”(빌리 조엘)에 따라 살면 결코 과거를 잊지 못할 것이며 “절대 말하지 않을 것들”(에이브릴 라빈)의 목록이 길수록 인생에서 원하는 것을 얻을 가능성이 줄어든다는 것을 조심하라고 말합니다. 따라서 다른 사람이 “나랑 얘기 좀 하자”(아니타 베이커)는 요청을 할 때는 경계심을 풀고 마음속에 있는 이야기를 “소통”(B-52)하는 것이 중요합니다. 하지만 적어도 일부 업무 관련 상황에서는 반드시 그렇지는 않다고 M. S. Cole 외(2008)의 팀워크에 관한 연구에 따르면 말합니다. 이 연구자들에 따르면, 때로는 함께 일하는 다른 사람들이 자신을 괴롭히는 행동이나 말에 대해 자신의 감정을 숨기는 것이 더 나을 수 있으며, 그러한 감정이 팀의 관심의 초점이 되어 팀이 적시에 효율적인 방식으로 작업을 수행하는 데 방해가 되지 않도록 하는 것이 더 나을 수 있습니다. 이 연구는 조건부 프로세스 모델의 추정과 해석의 메커니즘을 설명하는 첫 번째 사례의 데이터를 TEAMS라는 데이터 파일로 제공하며, www.afhayes.com 에서 확인할 수 있습니다. 이 연구는 자동차 부품 제조 회사에 고용된 60개의 작업 팀을 대상으로 진행되었으며, 회사 직원 200여 명을 대상으로 작업 팀에 대한 일련의 질문과 팀 감독자에 대한 다양한 인식에 대한 설문조사에 대한 응답을 기반으로 합니다. 이 연구의 일부 변수는 그룹 수준에서 측정된 것으로 같은 팀원들이 말한 내용을 종합하여 도출한 것입니다. 다행히 팀원들이 팀에 대한 질문에 응답하는 방식이 매우 유사하여 이러한 종류의 집계를 정당화할 수 있었습니다. 다른 변수는 순전히 팀 상사의 보고를 기반으로 합니다.

이 분석과 관련된 네 가지 변수를 측정했습니다. 팀원들이 다른 팀원들의 업무를 약화시키거나 변화와 혁신을 방해하는 행동을 얼마나 자주 했는지 등 팀원들의 역기능적 행동에 대한 일련의 질문(데이터 파일에서 DYSFUNC, 점수가 높을수록 팀 내 역기능적 행동이 많음을 나타냄)을 던져 팀원들의 역기능적 행동을 측정했습니다. 또한 팀원들에게 직장에서 ‘화가 났다’, ‘역겨웠다’ 등을 얼마나 자주 느끼는지 물어봄으로써 그룹의 부정적인 정서적 분위기를 측정했습니다(NEGTONE, 점수가 높을수록 업무 환경의 부정적인 정서적 분위기를 더 많이 반영함). 팀 상사에게는 팀이 얼마나 효율적이고 적시에 일을 처리하는지, 팀이 생산 목표를 달성하는지 등 전반적인 팀 성과에 대한 평가를 제공하도록 요청했습니다(데이터의 성과, 점수가 높을수록 성과가 좋음을 반영하는 척도). 또한 슈퍼바이저는 팀원들이 자신의 감정에 대해 보내는 비언어적 신호를 얼마나 쉽게 읽을 수 있는지를 측정하는 일련의 질문, 즉 비언어적 부정적 표현력(데이터 파일의 NEGEXP, 점수가 높을수록 팀원들이 부정적인 감정 상태를 비언어적으로 더 잘 표현한다는 의미)에 응답했습니다. 이 연구의 목표는 업무 팀원의 역기능적 행동이 업무 팀의 성과에 부정적인 영향을 미칠 수 있는 메커니즘을 조사하는 것이었습니다. 연구진은 역기능적 행동(X)으로 인해 상사와 다른 직원들이 직면하고 관리하려고 하는 부정적인 감정(M)으로 가득 찬 업무 환경이 조성되면 업무에 집중하지 못하고 업무 수행에 방해가 된다는 중재 모델을 제안했습니다(Y). 그러나 이 모델에 따르면 팀원들이 부정적인 감정(W)을 조절할 수 있게 되면, 즉 자신의 감정을 다른 사람에게 숨길 수 있게 되면 업무 환경의 부정적인 분위기와 다른 사람의 감정을 관리하는 데 집중할 필요 없이 당면한 업무에 집중할 수 있게 됩니다. 즉, 이 모델에서는 업무 환경의 부정적인 정서적 어조가 팀 성과에 미치는 영향은 팀원이 자신의 감정을 숨기는 능력에 따라 달라지며, 부정적 감정을 숨기지 않고 표현하는 팀에서 부정적인 정서적 어조가 성과에 미치는 부정적 영향이 더 강하다고 가정합니다.

Data: Introduction to Mediation, Moderation, and Conditional Process Analysis; “data files and code”

teams <- read_csv("data/hayes2022data/teams/teams.csv")

뒷 부분: \(Y \leftarrow X, M (W)\)

mod_m <- lm(perform ~ dysfunc + negtone + negexp + negtone:negexp, data = teams)
# 위와 동일: dysfunc은 통제 변수로 볼 수 있음
mod_m <- lm(perform ~ dysfunc + center(negtone) * center(negexp), data = teams)
S(mod_m, brief = T)
Coefficients:
                               Estimate Std. Error t value Pr(>|t|)   
(Intercept)                    -0.03207    0.05866  -0.547  0.58680   
dysfunc                         0.36606    0.17782   2.059  0.04429 * 
center(negtone)                -0.43144    0.13118  -3.289  0.00176 **
center(negexp)                 -0.04357    0.11343  -0.384  0.70239   
center(negtone):center(negexp) -0.51697    0.24092  -2.146  0.03632 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 0.4488 on 55 degrees of freedom
Multiple R-squared: 0.312
F-statistic: 6.235 on 4 and 55 DF,  p-value: 0.0003276 
  AIC   BIC 
80.92 93.49 

앞 부분: \(M \leftarrow X\)

mod_x <- lm(center(negtone) ~ dysfunc, data = teams)
S(mod_x, brief = T)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.02148    0.06176  -0.348 0.729182    
dysfunc      0.61975    0.16683   3.715 0.000459 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard deviation: 0.4763 on 58 degrees of freedom
Multiple R-squared: 0.1922
F-statistic:  13.8 on 1 and 58 DF,  p-value: 0.000459 
  AIC   BIC 
85.23 91.51 
predictorEffects(mod_m, ~negtone) |> plot()

avPlots(mod_m)

Hayes의 PROCESS를 매크로: Model 14

  • moments = 1: moderator의 세 값 (mean, mean +/- 1 SD)에 대한 조건적 효과를 계산
  • center = 1: centering
  • Index of moderated mediation: 조절변수가 1 변할 때, 간접효과의 변화가 얼마나 되는지를 나타내는 지수. 즉 간접효과가 얼마나 조절되는가?
process(data=teams, y="perform", m="negtone", w="negexp", x="dysfunc", jn=0, moments=1, center=1, model=14)

# ********************* PROCESS for R Version 4.3.1 ********************* 
 
#            Written by Andrew F. Hayes, Ph.D.  www.afhayes.com              
#    Documentation available in Hayes (2022). www.guilford.com/p/hayes3   
 
# *********************************************************************** 
               
# Model : 14     
#     Y : perform
#     X : dysfunc
#     M : negtone
#     W : negexp 

# Sample size: 60

# Random seed: 791081


# *********************************************************************** 
# Outcome Variable: negtone

# Model Summary: 
#           R      R-sq       MSE         F       df1       df2         p
#      0.4384    0.1922    0.2268   13.7999    1.0000   58.0000    0.0005

# Model: 
#              coeff        se         t         p      LLCI      ULCI
# constant   -0.0215    0.0618   -0.3479    0.7292   -0.1451    0.1021
# dysfunc     0.6198    0.1668    3.7148    0.0005    0.2858    0.9537

# *********************************************************************** 
# Outcome Variable: perform

# Model Summary: 
#           R      R-sq       MSE         F       df1       df2         p
#      0.5586    0.3120    0.2015    6.2350    4.0000   55.0000    0.0003

# Model: 
#              coeff        se         t         p      LLCI      ULCI
# constant   -0.0321    0.0587   -0.5467    0.5868   -0.1496    0.0855
# dysfunc     0.3661    0.1778    2.0585    0.0443    0.0097    0.7224
# negtone    -0.4314    0.1312   -3.2890    0.0018   -0.6943   -0.1686
# negexp     -0.0436    0.1134   -0.3841    0.7024   -0.2709    0.1838
# Int_1      -0.5170    0.2409   -2.1458    0.0363   -0.9998   -0.0341

# Product terms key:
# Int_1  :  negtone  x  negexp      

# Test(s) of highest order unconditional interaction(s):
#       R2-chng         F       df1       df2         p
# M*W    0.0576    4.6043    1.0000   55.0000    0.0363
# ----------
# Focal predictor: negtone (M)
#       Moderator: negexp (W)

# Conditional effects of the focal predictor at values of the moderator(s):
#      negexp    effect        se         t         p      LLCI      ULCI
#     -0.5437   -0.1504    0.2129   -0.7063    0.4830   -0.5770    0.2763
#      0.0000   -0.4314    0.1312   -3.2890    0.0018   -0.6943   -0.1686
#      0.5437   -0.7125    0.1530   -4.6565    0.0000   -1.0192   -0.4059

# *********************************************************************** 
# Bootstrapping progress:
#   |>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>| 100%

# **************** DIRECT AND INDIRECT EFFECTS OF X ON Y ****************

# Direct effect of X on Y:
#      effect        se         t         p      LLCI      ULCI
#      0.3661    0.1778    2.0585    0.0443    0.0097    0.7224

# Conditional indirect effects of X on Y:

# INDIRECT EFFECT:

# dysfunc    ->    negtone    ->    perform

#      negexp    Effect    BootSE  BootLLCI  BootULCI
#     -0.5437   -0.0932    0.1537   -0.3751    0.2613
#      0.0000   -0.2674    0.1193   -0.5282   -0.0591
#      0.5437   -0.4416    0.1610   -0.7788   -0.1448

#      Index of moderated mediation:
#            Index    BootSE  BootLLCI  BootULCI
# negexp   -0.3204    0.1887   -0.7753   -0.0491

# ******************** ANALYSIS NOTES AND ERRORS ************************ 

# Level of confidence for all confidence intervals in output: 95

# Number of bootstraps for percentile bootstrap confidence intervals: 5000

# W values in conditional tables are the mean and +/- SD from the mean.
 
# NOTE: The following variables were mean centered prior to analysis: 
#          negexp negtone

Exercises

A. 다음 데이터는 the Octogenarian Twin Study of Aging에서 나타나는 패턴을 기반으로 생성한 데이터

출처: Longitudinal Analysis: Modeling Within-Person Fluctuation and Change by Lesa Hoffman

includes 550 older adults age 80 to 97 years.
Cognition was assessed by the Information Test, a measure of general world knowledge (i.e., crystallized intelligence; range = 0 to 44)
demgroup 1: those who will not be diagnosed with dementia (none group = 1; 72.55%),
demgroup 2: those who will eventually be diagnosed with dementia later in the study (future group = 2; 19.82%)
demgroup 3: those who already have been diagnosed with dementia (current group = 3; 7.64%)

spss_chapter2.csv

  • 우선 demgroupsexMW를 factor로 변환하고,
  • 선형모형 cognition ~ grip * age + demgroup * sexMW을 세우고, 해석해보세요.
cognition <- read_csv("data/spss_chapter2.csv")
cognition |> head(5)
# A tibble: 5 × 6
  PersonID cognition   age  grip sexMW demgroup
     <dbl>     <dbl> <dbl> <dbl> <dbl>    <dbl>
1        1        23  92.6     9     1        1
2        2        24  91.8    11     0        2
3        3        29  92.6    12     0        1
4        4        16  94.4     6     1        1
5        5        27  85.8     9     1        1

B. 다음 데이터는 TV 시청시간(tv)과 개인적인 능력(ability)이 학생들의 학업성취도(achieve)에 미치는 영향을 알아보기 위한 연구의 일부입니다.

Source: Multiple Regression and Beyond (3e) by Timothy Z. Keith
Data: tv_ability.csv

  • TV 시청시간(tv)이 학업성취도(achieve)에 미치는 부정적인 영향이 개인적인 능력(ability)에 따라 차등적으로 발생한다는 가설을 테스트하기 위한 분석을 해보세요. 즉, 개인적인 능력이 moderator가 되는지를 알아보세요.
  • SES(ses) 변수를 통제할 필요성이 있는지 살펴보고, 통제한 경우와 하지 않은 경우를 비교해보세요.
tv <- read_csv("data/tv_ability.csv")
tv
# A tibble: 500 × 4
     ses ability    tv achieve
   <dbl>   <dbl> <dbl>   <dbl>
 1 -0.2       96     3      41
 2  1.15     107     4      55
 3  0.1       97     2      46
 4 -1.51      84     5      38
 5  1.87     116     1      56
 6  1.15     102     2      42
 7 -1.27      97     7      44
 8 -0.93      95     5      55
 9 -0.37      87     5      42
10 -1.98      85     6      40
# ℹ 490 more rows