왜 조건적 가설을 수립하고 검증하기 위해서는 상호작용 모델을 사용해야만 하는가?

이론적으로 가산적 관계를 상정한 선형 모델은 각각의 예측변수들이 서로 독립적으로 종속변수와 관계를 맺고 있을 것이라고 봅니다. 이는 선형 모델의 경우 각 예측변수에 대해 조건적 가설이 아닌 독립적 가설을 수립하고 있다고 이해할 수 있습니다. 논문에 보시면 흔히, “다른 조건들이 모두 일정할 때, \(x\)의 \(y\)에 대한 효과는 어떠할 것이다” 라는 형태의 진술을 가설로 사용하는 것이 그러합니다. 다른 조건들이 모두 일정하다는 얘기는 다중선형회귀모델에서 우리가 관심을 가지고 있는 \(x\) 이외의 다른 모든 변수들을 상수로 고정하였을 때, \(x\)가 \(y\)에 미치는 부분적 효과(partial effect)만을 보겠다는 것이고, 이것이 이 챕터 맨 처음에 제시하였던 식에서의 \(\beta\)라고 이해할 수 있습니다.

그러나 여기서는 정치학이라고도 할 수 있겠습니다만 사회과학 제분야의 경우 어떠한 맥락 혹은 조건의 효과를 고려하는 것이 모델링에 주요한 영향을 미치기 때문에 상호작용항을 포함한 모델을 생각해보아야 하는 경우가 생길 때가 많습니다. 예를 들어, 한국의 선거정치 속에 나타나는 계급투표에 대해 관심이 있다고 해보겠습니다. 이때, 기존 연구들의 주요 가설은 소득 수준에 따라 돈 많은 사람은 세금 많이 떼기 싫고 기득권 층일테니 보수정당을 지지할 것이라고 기대하고, 반대로 돈 없는 사람은 복지지출 등을 제고하기를 기대하며 진보정당에 투표할 가능성이 클 것이라고 해보겠습니다. 이때, 나의 연구가설은 이러한 계급적 투표가 실제 선거에서 정당 선택으로 이어지는 관계가 개별 유권자들의 정치적 세련도, 혹은 정치지식 수준에 따라 달라진다는 조건적 가설일 수 있습니다. 왜냐하면 돈이 없더라도 정치지식 수준이 높은 사람들의 경우 특정 정책에 대한 이해도가 더 높아 정당 선호가 달라질 수 있으니까요. 다양한 논의가 가능하겠지만 일단 이런 맥락에서 사회과학 분야에서 상호작용 모델은 상당히 빈번하게 논의됩니다.

마찬가지로 Brambor al. (2006)도 조건적 가설의 사례를 하나 제시하는데, 그들이 제시하는 가설은 조금 더 통계적이라고 해야할까요, 도식적입니다. 실제 연구문제를 바탕으로 한 가설은 아닙니다. 일단 조건적 가설이라고 보기 어려운 경우를 한 번 보겠습니다.

그들은 “\(Z\)라는 조건이 존재할 때, \(X\)의 한 단위 증가가 \(Y\)와 관계를 가지고, \(Z\)라는 조건이 부재할 경우에는 \(X\)와 \(Y\) 간의 관계 또한 성립되지 않을 것이다”라는 가설을 제시합니다. 즉, 여기서 \(Z\)라는 변수는 특정 조건의 존재 여부를 보여주므로 이항변수(binary variable)이라고 할 수 있겠죠? 만약 우리가 이항변수 \(Z\)를 성별이라고 한다면, 이 성별 변수의 계수값은 \(Y\)에 있어서 남성과 여성일 경우 각기 다른 절편값으로 해석할 수 있을 겁니다. 남성이 1, 여성이 0이라고 한다면

\[\begin{equation*} \begin{aligned} &\text{여성}: Y = b_0 + b_1X + b_2(Z=0) + \epsilon = b_0 + b_1X + \epsilon \\ &\text{남성}: Y = b_0 + b_1X + b_2(Z=1) + \epsilon = b_0 + b_1X + b_2 + \epsilon \end{aligned} \end{equation*}\]

이때 이 두 식의 차이는 \(b_2\), 상수로 나타납니다. 따라서 이 경우에는 조건적 가설이라고 부르기가 어렵습니다.

왜 상호작용 모델에 모든 구성항을 다 포함해야할까?

그렇다면 왜 상호작용 모델에 모든 구성항을 다 포함해야만 할까요? 그리고 상호작용 모델에서의 구성항과 일반 가산적 관계의 선형 모델에서의 예측변수 간의 차이는 무엇일까요?

상호작용 모델(\(Y = b_0 + b_1X + b_2Z + b_3XZ + \epsilon\))에서 구성항을 제외한다면, 우리는 \(Y = b_0 + b_1XZ + \epsilon\)라는 모델을 생각해볼 수 있습니다. 구성항을 제외해도 말이 된다라는 얘기는 위의 \(Y = b_0 + b_1X + b_2Z + b_3XZ + \epsilon\)이 \(Y = b_0 + b_1XZ + \epsilon\)의 모델과 다르지 않다라는 주장을 해야한다는 것을 의미합니다. 이때, 만약 구성항을 제외한 모델이 타당하다면,

1. \(Z\)가 평균적으로 \(Y\)에 미치는 효과가 존재하지 않다거나

2. \(X\)가 0일 때, \(Z\)가 \(Y\)에 미치는 효과가 없다

라고 주장할 수 있어야만 합니다. 다시 말하면, 첫째로 \(Z\)가 평균적으로 \(Y\)에 미치는 효과가 없다는 얘기는 사실상 비체계적 요인으로서 \(Y\)에 평균적으로 체계적 변화를 가져오는 요인이 \(X\) 뿐이라고 주장하던가(이 경우 본질적으로 해당 모형의 함의는 \(Y = b_0 + b_1X + \epsilon\)과 다를 바 없게 됨), 혹은 \(X\)가 0일 때, \(b_1XZ\) 항이 제거되면서 절편값인 \(b_0 + \epsilon\)만 남게 되어 두 번째 주장이 타당하다고 말할 수 있어야 합니다. 그러나 위의 두 주장은 가설로서 거의 정당화되기 어렵습니다(Brambor et al. 2006, 66). 왜냐하면 \(Y = b_0 + b_1X + b_2Z + b_3XZ + \epsilon\)에서 \(b_2\)는 상호작용 모델에서 \(X\)가 0일때의 \(Z\)가 \(Y\)에 미치는 효과를 보여주기 때문에, 이는 \(Z\)가 평균적으로 \(Y\)에 미치는 효과가 존재않다는 것을 보여주어야 하는 것과 배치됩니다. 마찬가지로 \(Z\)의 \(Y\)에 대한 평균적인 효과가 0이라고 하더라도 \(b_2\)가 반드시 0이라는 보장도 없습니다. \(X\)가 0일 때, \(Z\)가 \(Y\)에 대해 가지는 효과가 전혀 없다(\(b_2 = 0\))고 미리 전제하는 것보다는, \(b_2\)가 0이 아닐 수도 있을 가능성을 먼저 생각해보는 것이 더 많은 경우의 수를 제공합니다(Brambor et al. 2005, 66). 다시 말해, 선험적으로 구성항을 제외하는 분석은 \(b_2\)가 0이 아니라면 우리가 관심을 가지고 있는 모수의 추정치를 왜곡하여 잘못된 추론을 도출하게 할 수 있습니다(Brambor et al. 2006, 68).

위에서 언급한 바와 같이, \(X\)에 대한 계수, \(b_1\)와 \(Z\)에 대한 계수 \(b_2\)는 일반적인 가산적 관계의 모델(\(Y =b_0 + b_1X + b_2Z + \epsilon\))의 \(X\)와 \(Z\)의 계수들과는 다릅니다. 가산 모델에서 \(b_2\)는 \(Z\)의 \(Y\)에 대한 평균적인 효과를 보여줍니다. 그러나 상호작용 모델에서 \(b_2\)는 \(X\)에 따라 조건적으로 변화하는 \(Z\)의 \(Y\)에 대한 효과를 보여줍니다.

왜 실질적으로 유의미한 한계효과와 표준오차를 보여주어야 하는가?

상호작용 모델에 있어서 각 예측변수들이 종속변수에 미치는 효과를 제대로 보여주기 위해서 Brambor et al. (2006)은 실질적으로 유의미한 한계효과와 표준오차를 함께 제시해야 한다고 주장하고 있습니다. 그 이유는 상호작용항의 효과가 일정하지 않기 때문(not constant)입니다. 상호작용 모델(\(Y = b_0 + b_1X + b_2Z + b_3XZ + \epsilon\))에서 \(X\)의 한 단위 증가와 관계된 \(Y\)의 변화분을 살펴보기 위해 편미분을 할 경우, 우리는 아래와 같은 식을 얻게 됩니다.

\[ \frac{\partial Y}{\partial X} = b_1+ b_3 Z \]

위의 식은 상호작용 모델에서 \(X\)의 \(Y\)에 대한 효과가 \(Z\) 변수의 값에 따라서 조건적으로 변화한다는 것을 의미합니다 (Brambor et al. 2006, 73). 따라서 상호작용항의 계수는 직접적으로 해석될 수는 없습니다. 그 값이 \(Z\)에 따라 변화하니까요. 우리는 계수의 부호나 통계적 유의성에 대해서는 이야기할 수 있지만 실질적으로 그 효과의 크기 등에 대해서는 계수만 보고는 대답할 수 없게 됩니다.

때문에 \(Z\)에 의해 조건적으로 변화하는 \(Y\)에 대한 \(X\)의 효과를 살펴보기 위해서는 \(X\)의 \(Y\)에 대한 한계효과를 계산할 필요가 있습니다. 한계효과는 \(Z\)의 조건 하에서 \(X\)의 한 단위 변화가 평균적으로 얼마만큼의 \(Y\)의 변화와 관계되는지를 보여줍니다. 또한 한계효과의 표준오차는 그렇게 계산된 한계효과가 얼마나 확실한지를 보여줍니다. 한계효과의 표준오차가 가지는 함의는 일반적으로 선형회귀모델의 계수값과 표준오차 간의 관계와 비슷하다고 이해하셔도 무방합니다. 그 한계효과가 통계적으로 유의미한 효과인지를 보여주는 것이니까요.

정리하자면 Brambor et al. (2006)의 함의는 다음과 같습니다.

• 상호작용 모델이라고 하더라도 상호작용항뿐만 아니라 모든 구성항을 포함시켜 분석해야 한다.

• “\(Z\)를 상수로 고정한다(=통제한다)”는 것은 \(Z\)가 0이라는 것과 같은 의미가 아니다.

• 모델 내에 존재하는 두 구성항의 곱으로 상호작용항을 만들었기 때문에 다중공선성(multicollinearity)이 발생할 수 있다.²

• 상호작용항의 계수는 \(\frac{\partial Y}{\partial X} = b_1+ b_3 Z\)로 나타낼 수 있고, 이때 상호작용의 계수가 편미분을 하더라도 \(Z\)라는 새로운 변수에 의해 조건적으로 변화할 수 있다는 것을 알 수 있다. 따라서 상호작용항의 계수를 직접적으로 일반선형회귀모델의 계수처럼 해석하기는 어렵다.

그렇다면 한 번 실제 데이터를 통해 분석해보도록 하겠습니다.

기본 아이디어

우리가 알고싶은 것은 모집단의 특성입니다. 그러나 모집단을 관측하거나 혹은 모집단에 관한 완전한 데이터를 얻을 수 없기 때문에 우리는 모집단의 특성을 잘 반영할 것이라고 기대된 표본을 통해 모집단이 어떠할 것이라는 통계적인 추론을 하게 됩니다.

• 즉, 추론의 기저에는 표본이 모집단을 잘 대표하고 있다면 표본에서 우리가 관측한 분포, 특성들이 모집단에서도 그러할 것이라고 주장하는 것이 가능하다는 가정이 깔려 있습니다.

이론적으로 우리는 모집단에서 무수히 많은 표본을 추출할 수 있고, 우리는 이렇게 반복된 표본들로 연구를 반복할 수 있습니다.

• 그러나 모집단으로부터 표본을 추출하는 표집방법의 근본적인 문제로 동일한 모집단에서 뽑은 표본들이라고 할지라도 표본들 간에는 약간의 차이가 존재할 수 있습니다.

선형회귀모델을 한 번 생각해보겠습니다. 우리는 기본적으로 함수적 관계를 통해 어떠한 예측변수들이 종속변수와 관계를 맺을 것이라고 주장합니다. 그러나 약간씩 다른 표본의 자료들을 모델에 투입하는 순간, 당연히 결과로 추정되는 계수값들도 차이가 나게 될 것입니다. 우리는 이렇게 서로 다른 표본들로부터 얻은 통계치, 계수값들을 한 데 모아 그 분포를 그려볼 수 있습니다. 이것이 바로 표집분포입니다.

가우스-마르코프 가정에 따르면, 우리는 이렇게 한 데 그러모은 표본들로부터 나온 계수들의 기대값이 모집단 수준의 모수, 단 하나의 진실된 값과 동일할 것이라고 기대합니다.

• 그러나 현실적으로 모집단으로부터 수없이 많은 표본들을 얻는 것은 불가능합니다.

• 예를 들어, 정치학에서 흔히 사용되는 데이터셋인 POLITY IV를 생각해보겠습니다.

  • POLITY IV는 어디까지나 측정가능한 민주주의 정치체제에 대한 제반 속성들을 수치화한 것입니다. 민주주의라고 하는 모집단을 의미하는 것은 아니죠. 어디까지나 표본입니다.

  • 하지만 이 데이터셋 하나를 구축하는 데만도 엄청난 자원과 비용이 소요되었습니다. 네, 현실에서는 이렇게 표본 하나를 제대로 구하는 것도 쉽지 않습니다. 그러므로 사실 우리는 단 하나의 표본에서 나타나는 특성과 통계치를 이용해서 표집분포를 추정하고 있습니다.

  • 이렇게 하기 위해서는 몇 가지 가정들이 충족되어야만 합니다. 예를 들어, 우리가 선형회귀모델을 단 하나의 표본만을 가지고 추정한다고 할 때, 오차항의 독립성과 정규분포가 가정되어야 합니다. 만약 그렇지 않다면 우리는 추정한 계수값이 일관되거나 편향되지 않았다고 하기 어려울 것입니다.

부트스트랩 방법은 표집분포를 추정하는 데 있어서 다른 접근법을 취합니다. 부트스트랩 방법은 대개 복원추출을 통해 모집단으로부터 얻은 \(n\)개의 관측치를 가진 표본으로부터 다시 \(n\)개의 관측치를 가진 표본을 뽑아냅니다.

• 즉, 표본에서 표본을 뽑습니다. 이렇게 재추출된 표본들은 원래 모집단에서 추출한 표본과 동일한 값을 지닐수도, 아닐수도 있습니다. 왜냐하면 무작위 복원추출을 했기 때문에 어떤 값은 중복되고 어떤 값은 재추출된 표본에는 포함되지 않을수도 있으니까요.

• 부트스트랩은 이 과정을 반복합니다. 그러고나면 우리는 모집단에서 얻은 단 하나의 표본으로부터 무수히 많은 재추출 표본들을 얻을 수 있게 됩니다.

왜 이런 방법을 사용할까요? 목적은 다르지 않습니다. 우리는 여전히 우리가 알 수 없고, 관측할 수 없는 모집단의 특성을 알고 싶습니다. 그 모집단의 특성을 알기 위해서는 그 특성의 확률분포를 알아야할텐데, 우리는 그 분포가 어떠할 것이라는 확신을 가질 수가 없습니다.

• 즉, 우리는 사전에 모집단의 확률분포가 어떠할 것이라고 알지 못합니다. 우리는 오직 관측된 표본 통계치들만을 가지고 모수를 추정할 뿐입니다. 한정된 데이터로 인한 불확실성은 우리의 추론을 제약합니다.

• 또한, 대부분의 통계방법들은 표본의 크기에 따라서 다른 결과들을 내어놓기도 합니다. 표본의 크기 자체가 커질수록 우리는 모수에 대한 더 안정적인 결과를 얻을 수 있을 것이라 기대합니다. 하지만 비용과 시간의 제약으로 인하여 표본 규모를 늘리기란 쉽지 않습니다. 작은 규모의 표본을 가지고 할 수 있는 대안이 없을까? 여기서 부트스트랩 방법이 시작된 것입니다.

• 부트스트랩 방법은 기본적으로 하나의 표본으로부터 수많은 부트스트랩 표본들을 추출함으로서 앞서의 문제들을 극복하고자 하며, 나아가 하나의 표본 특성을 부트스트랩 표본들로부터 얻은 표집분포로 추론하고, 나아가 그것을 모집단의 모수를 추론하는데까지 사용하고자 합니다.

부트스트랩 방법을 사용하게 되면 우리는 무수히 많은 부스트스랩 표본들을 하나의 표본으로부터 얻게 됩니다. 그리고 부트스트랩 표본들은 하나의 표본으로부터 무작위 복원추출을 한 결과물입니다. 이제 기본적인 논리가 이해가 가실 겁니다.

• 부트스트랩 표본들 \(\rightarrow\) 표본 \(\rightarrow\) 모집단으로 가는 경로인 것이죠.

• 부트스트랩 방법은 원 표본(original sample)과 동일한 규모의 표본으로 재추출하고, 부트스트랩 표본들로부터 얻은 표본 통계치의 분포를 표집분포로 사용합니다.

  • 즉, 부트스트랩 방법은 원 표본을 일종의 모집단에 대한 대리(proxy)로 간주하고 무작위 반복 표집을 하는 것이죠.

  • 이렇게 새로 뽑은 부트스트랩 표본은 흡사 원표본에 대한 ’하위 표본’으로 보이지만 정확하게는 ’하위 표본’은 아닙니다. 모집단에 대한 또 다른 표본이죠. 이것이 가능하기 위해서는 하나의 가정이 필요합니다.

  • 바로 원 표본이 모집단을 대표할 수 있는 무작위 반복추출된 표본이어야 합니다. 어디까지나 원 표본이 잘 추출되었다면, 거기서 무작위로 반복추출해서 다시 만든 부트스트랩 표본들은 모집단에서 무작위 반복추출한 결과와 다르지 않을 것이라는 굉장히 강력한 무작위화(randomization)에 대한 믿음이 뒷받침하고 있습니다.

• 결과적으로 부트스트랩의 재표집 과정은 무수히 많은 표본들을 얻을 수 있고, 부트스트랩 표본들로부터 얻은 표본통계치들은 동일한 모집단에서 추출된 무작위 표본들 간의 차이(variability)를 추정하는 데 사용될 수 있습니다.

이러한 부트스트랩 표집분포를 이용해 우리는 실제 관측치를 이용한 신뢰구간 측정이나 가설검정을 수행할 수 있게 됩니다.

• 여기서 한 가지 착각하면 안 되는 것은 비모수 부트스트랩 방법이 결코 모집단의 분포에 대한 어떠한 가정을 하고 있는 것은 아니라는 점입니다.

• 우리는 작은 규모의 표본을 가지고 있을 때, 이 표본이 최소자승법을 사용하기 위한 가정들을 충족시키는지 여부를 알지 못하기 때문에 OLS의 추정 결과를 얻더라도 그것이 과연 BLUE인지를 검정할 방법이 없습니다.

• 하지만 부트스트랩 방법을 이용하면 이 문제를 시뮬레이션한 표집 분포를 이용해 해결할 수 있습니다. 이것이 소규모 표본을 이용할 때 표본 통계치의 신뢰도⁴가 낮아 생길 수 있는 문제들을 다룰 수 있는 부트스트랩의 장점입니다.

부트스트랩 방법의 장점은 우리가 작은 표본을 가지고 있거나 모집단의 분포를 모른다고 하더라도 모집단의 특성을 추정할 수 있도록 해준다는 것입니다.

• 시뮬레이션된 표집분포를 가지고 우리는 계수, 표준오차, 그리고 신뢰구간 등 기존에 추정하던 통계치들을 모두 추정할 수 있습니다. 또한 OLS가 회귀선을 그리는 방법을 부트스트랩 방법과 비교하자면, 부트스트랩 방법은 사전에 모집단의 분포가 어떠할 것이라고 가정을 하지 않기 때문에 모집단의 확률분포가 정규분포가 아니라고 가정하는 다른 여러 추정기법들과도 함께 사용될 수 있습니다.

• 개인적으로는 부트스트랩 방법은 대개 소규모 표본의 문제로 여러가지 어려움을 겪는 사회과학 분야의 경험연구에 있어서 특히나 유용한 대안이라고 생각합니다. 그러나 분명 부트스트랩 방법에도 한계는 있습니다.

  • 바로 무작위화의 가정이 깨어질 경우입니다. 만약 우리가 가진 단 하나의 표본이 모집단으로부터 무작위 복원추출되었다고 할 근거가 충분하지 못하다면 어떻게 될까요? 원 표본이 편향된 순간, 부트스트랩 표본들도 모집단의 특성을 대표하는 것이 아니라 사실상 원 표본의 체계적 편향을 반영한 표본들이 되고 맙니다.

부트스트랩… 어디다 써먹지?

정리하자면 비모수 부트스트랩을 통해 우리는 표집분포를 얻을 수 있습니다. 그리고 이 표집분포에 대한 함수적 관계를 상정하는 것도 자유롭습니다.

• 분석적 접근법(analytic approach)에 따르면 모델 \(y = \beta_0 + \beta_1x + \beta2_z + \beta_3xz + \epsilon\)에서 상호작용을 나타내는 편미분 결과 \(\hat{\beta_1} + \hat{\beta_3}z\)의 표준오차, 즉 표본의 차이에 따라 그 상호작용 효과가 다르게 나타날 수 있는 변동성은 \(Var(\hat{\beta_1}) + z^2Var(\hat{\beta_3}) + 2z Cov(\hat{\beta_1}, \hat{\beta_3})\) 으로 나타납니다.

• 하지만 만약 우리가 \(g = 1, \cdots, G\)개의 부트스트랩 표본들을 추출할 수 있게 된다면, 그 부트스트랩 표본들로부터 얻은 \(\hat{\beta}_1^{(1)} + \hat{\beta}_3^{(1)}z, \cdots, \hat{\beta}_1^{(G)} + \hat{\beta}_3^{(G)}z\)를 가지고 요약통계를 보여주면 됩니다. 훨씬 간단하죠.

행렬 OLS를 통해 우리가 얻을 수 있는 예측값이 \(\hat{\textbf{y}} = \bar{\textbf{X}}\hat{\beta}\)이라는 것을 다시 떠올려보겠습니다.

• 부트스트랩 표본을 \(g=1,\cdots,G\)개 가졌다고 한다면 이 부트스트랩 표본들로부터 얻은 계수값은 \(K\times1\)의 벡터로 나타낼 수 있습니다: \(\hat{\beta}^{(g)}\). 따라서 하나의 부트스트랩 표본에서 얻어진 예측값은 \(\hat{\textbf{y}}^{(g)} = \bar{\textbf{X}}\hat{\beta}^{(g)}\)라고 할 수 있습니다.

• 그리고 \(G\)개의 부트스트랩 표본들을 이용해 알고 싶은 예측값에 대한 전체 부트스트랩 분포를 구할 수 있습니다. 간단히 말하면, \(x_1\) 변수에 대한 \(\hat{beta_1}\)은 부트스트랩 표본들 전체를 모델에 집어넣고 나온 \(\hat{\beta}^{(g)}\)의 분포로 나타낼 수 있다는 것입니다.

• 부트스트랩 표본들로부터 추정된 일련의 계수값들을 \(\hat{\textbf{b}} = (\hat{\beta}^{(1)} \cdots \hat{\beta}^{(G)})\)라고 하겠습니다.

• 그리고 부트스트랩 표본들에 모델을 적용하여 얻은 일련의 예측값들을 \(\hat{\textbf{Y}} = (\hat{y}^{(1)} \cdots \hat{y}^{(G)})\)라고 해보겠습니다.

• \(\hat{\textbf{Y}} = \bar{\textbf{X}}\hat{\textbf{b}}\)는 우리가 관심을 가지는 특정한 \(M\)이라는 예측값 각각에 대한 \(G\)개의 부트스트랩 표본들이라고 할 수 있습니다.

부트스트랩 R로 보기: 비모수 부트스트랩의 추출횟수와 결과

실제 데이터를 가지고 부트스트랩 방법을 적용해보겠습니다. 논리대로라면 부트스트랩 표본의 추출 횟수가 증가할수록, 부트스트랩 표본을 통해 얻은 계수 추정치들은 우리가 표집분포를 이용해 도출한 분석적 접근법의 결과에 근사하게 될 것입니다.

Boot <- QOG.s %>% 
  dplyr::select(wdi_gdpcapcon2010, p_polity2, wdi_trade)

# 비모수 부트스트랩 표본추출 횟수 변화에 따른 결과의 차이를 보겠습니다.
# 먼저 100번 부트스트랩을 할 경우입니다.
holder100 <- matrix(NA, 100) ## 일단 100번 추출한 결과가 들어갈 빈 행렬을 만듭니다.
for(i in 1:100){  ## 함수를 만듭니다. i는 1부터 100까지를 의미합니다.
  ind <- sample(1:nrow(Boot), ## ind는 전체 관측치(각 행의 번호) 중에서 
              ## QOG.s와 동일한 표본규모로 표집됩니다.
                size=nrow(Boot), replace=T) ## 복원표집된 결과입니다.
  # 이 ind가 헷갈리실 텐데, 뭐냐면 넘버링입니다. 매번 무작위로 넘버링이 
  # 나올거고 이 넘버링에 속하는 관측치들은 개별 부트스트랩 표본을 구성합니다.
  # 그 과정이 i번, 즉 1~100까지 반복되어 총 100개의 부트스트랩 표본이 생깁니다.
  mod <- lm(log2(wdi_gdpcapcon2010) ~ ## 모델은 앞서와 동일하게 적용합니다.
              1 +  p_polity2 + wdi_trade, 
            data=Boot[ind,]) ## ind에 속하는 행의 관측치만을 가지고 모델분석
  holder100[i,] <- coef(mod)[2] ## 이렇게 새로 부트스트랩 표본으로 모델을 
  # 분석해 나온 계수값만을 아까 만든 깡통행렬에 차곡차곡 하나씩 저장.
  # 민주주의에 대한 계수의 결과만을 보겠습니다.
}
head(holder100) ## 그려면 이렇게 총 100개의 계수값이 담긴 행렬이 완성됩니다.

##            [,1]
## [1,] 0.16737332
## [2,] 0.11392406
## [3,] 0.05170454
## [4,] 0.09482248
## [5,] 0.05588331
## [6,] 0.08455292

이 다음에는 단순하게 추출 횟수만 바꾸어주면 되겠죠? 빠르게 가겠습니다.

## 1000번 Bootstrapping
holder1000 <- matrix(NA, 1000)
for(i in 1:1000){
  ind <- sample(1:nrow(Boot), 
                size=nrow(Boot), replace=T)
  mod <- lm(log2(wdi_gdpcapcon2010) ~ 
              1 +  p_polity2 + wdi_trade,
            data=Boot[ind,])
  holder1000[i,] <- coef(mod)[2]
}

## 10000번 Bootstrapping
holder10000 <- matrix(NA, 10000)
for(i in 1:10000){
  ind <- sample(1:nrow(Boot),
                size=nrow(Boot), replace=T)
  mod <- lm(log2(wdi_gdpcapcon2010) ~ 
              1 +  p_polity2 + wdi_trade,
            data=Boot[ind,])
  holder10000[i,] <- coef(mod)[2]
}

## 100000번 Bootstrapping
holder100000 <- matrix(NA, 100000)
for(i in 1:100000){
  ind <- sample(1:nrow(Boot), 
                size=nrow(Boot), replace=T)
  mod <- lm(log2(wdi_gdpcapcon2010) ~ 
              1 +  p_polity2 + wdi_trade, 
            data=Boot[ind,])
  holder100000[i,] <- coef(mod)[2]
}

이렇게 계산된 각각의 100번, 1,000번, 10,000번, 100,000번의 부트스트래핑 결과로 얻어진 계수값들을 하나의 데이터로 만들어서 비교해보도록 하겠습니다. tidyverse 패키지의 tibble을 이용하여 하나의 데이터로 만들겠습니다.

# 각각의 부트스트랩한 계수값 결과를 하나의 데이터로 만들어줍니다.
npbs <- bind_rows(
  holder100 %>%
    as_tibble() %>% mutate(NPBS = "100") %>% 
    rename(Estimates = V1),
  holder1000 %>%
    as_tibble() %>% mutate(NPBS = "1000") %>% 
    rename(Estimates = V1),
  holder10000 %>%
    as_tibble() %>% mutate(NPBS = "10000") %>% 
    rename(Estimates = V1),
  holder100000 %>%           
    as_tibble() %>% mutate(NPBS = "100000") %>% 
    rename(Estimates = V1)) %>%
  mutate_if(is.numeric, round, 3)

# 각 부트스트랩 횟수별로 평균을 구해주겠습니다. 그리고 비교대상인 OLS의
# 민주주의 계수값 결과도 추가해주고요.
analytic <- lm(log2(wdi_gdpcapcon2010) ~ 
              1 +  p_polity2 + wdi_trade, data = QOG.s)

npbs <- npbs %>% group_by(NPBS) %>% mutate(
  `NPBS Mean` = mean(Estimates, na.rm = T),
  `OLS estimates` = coef(summary(analytic))[2,1],
  `OLS se` = coef(summary(analytic))[2,2]
)

부트스트랩으로 얻은 회귀계수와 OLS로 추정한 회귀계수를 비교해보도록 하겠습니다.

• 각 패널은 부트스트랩 복원표본추출 횟수를 보여줍니다.

• 그리고 x-축은 회귀계수 값을, y-축은 밀도입니다. 좌측 상단의 패널은 부트스트랩 패널을 100번 추출했을 때의 결과입니다.

• 붉은 선이 부트스트랩으로 얻은 회귀계수의 평균이고 소수점 셋째 자리에서 0.084입니다.

• 그리고 붉은 색으로 표시된 분포와 히스토그램이 부트스트랩으로 얻은 계수 관측치들의 분포를 보여줍니다.

• 검정색 선으로 표시된 분포는 OLS 분석으로 얻은 회귀계수값을 평균으로 하고 표준오차를 표준편차로 하는 정규분포입니다. 말이 되죠? 왜냐면 OLS의 가우스-마르코프 가정이 충족된다면 OLS의 회귀계수는 BLUE일테니, 그렇게 얻은 회귀계수의 기대값인 평균은 모집단의 모수, \(\beta\)와 같을 것이고 표집분포에서의 표준편차는 표본추출로 인한 표본들 간에 나타날 수 있는 차이를 보여주는 표준오차로 나타나니까요.

일단 부트스트랩 표본추출 100회차를 보면 벌써 평균이 0.002밖에 차이가 안 납니다. 다만 분포 양상은 정규분포라고 하기 조금 힘듭니다. 아무래도 관측치가 100개밖에 없다보니 좀 삐뚤빼뚤합니다. 히스토그램만 봐도 그렇죠.

두 번째로는 부트스트랩을 1,000번 시뮬레이션한 결과입니다. 평균 차이는 뭐 여전히 조금 존재하기는 하지만 분포가 한층 더 정규분포에 비슷해진 것을 확인할 수 있습니다.

• 마찬가지로 나머지 하단의 두 패널들도 점점 분포 양상이 정규분포 모양으로 수렴하는 것을 확인할 수 있습니다.

• 즉, 이를 통해서 우리는 부트스트랩 표본 추출 횟수가 증가할수록 부트스트랩 표본들을 이용해 얻은 계수값들의 분포는 OLS 회귀계수의 정규분포에 근사한다는 것을 알 수 있습니다. 시뮬레이션이 분석적 접근법의 대안일 수 있다는 것이죠.

부트스트랩 R로 보기(2): 신뢰구간과 점추정치

이전에 “변수에 대한 심화”라는 섹션에서 Gelman (2008)⁵에 따라 변수를 표준화하되 \(b=1sd(x)\)가 아니라 \(b=2sd(x)\)로 표준화를 해준 적이 있습니다.

• 결과적으로 \(b=2sd(x)\)로 표준화했을 때, 표준화 이전 계수의 경향성과 표준화한 이후의 계수의 경향성이 비슷한 양상을 보인다는 것을 확인했습니다.

• 어차피 표준화의 결과로 실질적 변수의 효과가 다르게 나타나는 것은 아니니만큼 \(b=2sd(x)\) 표준화가 직관적인 이해까지 도울 수 있다는 제언을 한 적이 있습니다.

이번에는 당시의 분석을 비모수 부트스트랩을 이용해서 재현해보도록 하겠습니다.

• 기억을 되살려보면 \(y_i = \delta_0 + \delta_1x_i, \:\:i\leq n\)이라고 할 때, \(a = \frac{1}{n}\sum_i x_i\)이고 \(b=2sd(x)\)인 \(y_i = \lambda_0 + \lambda_1\frac{x_i-a}{b}\)를 추정하는 것이었습니다.

• \(x\)에 맞추어 재구성하면 결국 아래와 같은 모델이라는 것을 알 수 있죠.

\[ y = \lambda_0 + \lambda_1(\frac{x_i - \bar{x}}{2\times\sigma_x}) \]

이번에는 QoG 데이터셋의 1인당 GDP와 민주주의 수준을 이용해 분석해보겠습니다. 데이터를 고려해보면 제 모델은 아래와 같습니다.

\[ \text{경제규모}_i = \lambda_0 + \lambda_1\frac{(\text{민주주의 수준})-(\text{민주주의 수준의 평균})}{\text{민주주의 수준의 표준편차}\times2} \]

# 2sd 표준화
QOG.std <- QOG.s %>%
  mutate(
    std.polity = p_polity2 - mean(p_polity2)/2*sd(p_polity2))

# 새롭게 모델링
model.sd <- lm(log2(wdi_gdpcapcon2010) ~ std.polity, data=QOG.std)

# 결과
model.sd %>% broom::tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% knitr::kable()

# 비교를 위한 원 변수 모델
model2 <- lm(log2(wdi_gdpcapcon2010) ~ p_polity2, data=QOG.std)
model2 %>% broom::tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% knitr::kable()

표준화를 했기 때문에 이제 예측변수는 측정단위에 따라 결과가 좌우되지는 않습니다. 즉, 변수가 여러개 있다면 표준화한 변수들의 효과들 간에는 상대적 비교가 가능합니다. 예를 들어, \(x_1\)이 \(x_2\)보다 \(y\)에 미치는 효과가 상대적으로 크다, 작다 등으로 해석이 가능하다는 얘기지요. 그것이 실질적인 효과와 직결되느냐는 조금 다른 문제입니다만 여하튼 표준화는 서로 분석단위가 \(kg\) 대 \(ml\) 등으로 상이한 변수들 간의 효과를 비교할 수 있도록 스케일링해주는 것입니다.

물론 여기서는 단 하나의 예측변수만을 사용하기 때문에 결과는 표준화된 예측변수—민주주의의 수준과 로그값을 취한 종속변수의 관계만을 보여줄 것입니다. 이 분석을 위하여 4,000번의 부트스트랩 복원표본추출을 해보도록 하겠습니다.

holder.std <- matrix(NA, 4000)
for(i in 1:4000){
  ind <- sample(1:nrow(QOG.std), size=nrow(QOG.std), replace=T)
  mod <- lm(log2(wdi_gdpcapcon2010) ~ std.polity, data=QOG.std[ind,])
  holder.std[i,] <- coef(mod)[2]
}

# 어떤 방식으로 데이터를 구성했는지 조금 직관적으로 이해할 수 있게
# 단순하고 늘어지는 코드를 써보았습니다.
# 먼저 부트스트랩 결과입니다.
std.data1 <- data.frame(ID="Bootstrapped",
                        Low=quantile(holder.std, 0.025), 
                        M=mean(holder.std),
                        High=quantile(holder.std, 0.975))
# 부트스트랩은 이미 우리가 4000개의 관측치를 가지고 있기 때문에 이 관측치에서
# 왼쪽꼬리 2.5백분위와 우측꼬리 2.5백분위에 해당하는 값이 95% 신뢰구간의
# 기준값이 됩니다.

# OLS 결과입니다.
std.data2 <- data.frame(ID="OLS",
                        Low=coef(summary(model.sd))[2,1] -
                          1.96*coef(summary(model.sd))[2,2],
                        M=coef(summary(model.sd))[2,1],
                        High=coef(summary(model.sd))[2,1] +
                          1.96*coef(summary(model.sd))[2,2])
std.data <- rbind(std.data1, std.data2)

자, 어떤가요? 부트스트랩 표본추출로 구한 계수 추정치와 OLS 추정치의 95% 신뢰구간 결과입니다. 보시다시피 부트스트랩 결과가 조금 큰 계수값과 넓은 신뢰구간을 가지고 있지만 두 계수 추정치는 매우 근사하고 이 차이는 실질적인 분석을 수행하는 데 있어서 무시할만합니다. 뭐 이렇게 말하면 조금 위험할수도 있지만요.

• 하지만 부트스트랩의 평균과 표준편차가 OLS의 분석적 결과인 계수 추정치와 표준오차 계산을 통한 신뢰구간보다 조금 과대추정되었더라도 부트스트랩 결과는 그 나름의 장점이 있습니다.

• 우리는 실질적으로 계수값들의 관측치를 통해 분포를 가지고 있기 때문에 그 분포를 이용해 원하는 값을 보여줄 수 있기 때문입니다. 게다가 부트스트래핑은 표본 규모가 작더라도 시도할 수 있기 때문에 어떤 부분에서는 OLS에 대해 상대적 이점이 있다고도 할 수 있습니다.

그렇다면 이번에는 비모수 부트스트랩을 넘어 이제 모수적 부트스트랩(parametric bootstrap), 흔히 정치학에서는 AJPS 2000년의 Gary King과 Michael Tomz, 그리고 Jason Wittenberg 논문으로 유명한 Clarify에 대해 살펴보도록 하겠습니다.

평균과 불확실성을 추정하기

실제 모델을 통해서 한 번 생각해보도록 하겠습니다. 이 모델은 임금수준과 교육수준, 그리고 그 두 변수의 상호작용을 통해 범죄율을 설명할 수 있다는 내용을 담고 있습니다. 정확히는 임금수준이 범죄율에 미치는 효과가 교육수준에 따라 조건적으로 나타날 것이라고 기대한다고 합시다.

\[ \text{범죄율} = \beta_0 + \beta_1\text{임금수준} + \beta_2{교육수준} + \beta_3(임금수\times교육수준) \]

여기서 우리가 관심있는 것은 바로 관찰된 교육 수준의 값이 평균일 때, 임금 수준의 변화에 따른 기대값이라고 할 수 있습니다.

• 모든 교육수준 관측치에 대해 임금 수준을 고임금(\(\overline{\text{임금}}\))과 저임금(\(\underline{\text{임금}}\))의 저임금으로 벡터화해보겠습니다.

• \(g = 1, \cdots. G\)라고 할 때, 다음과 같이 단계대로 진행합니다.

  1. \(\tilde{\beta}^{(g)} ~ N(\hat{\beta}, \hat{\sigma}^2(X'X)^{-1})\)의 분포를 따르는 \(\tilde{\beta}^{(g)}\)를 추출합니다.

  2. 교육수준 변수로부터 \(PE^{(g)}\)를 무작위로 추출합니다.

  3. \(\mu \overline{\text{임금}}^{(g)} = \tilde{\beta_0}^{(g)} + \tilde{\beta_1}^{(g)}\overline{\text{임금}} + \tilde{\beta_2}^{(g)}\text{교육수준}^{(g)} + \tilde{\beta_3}^{(g)}(\overline{\text{임금}}\times\text{교육수준}^{(g)})\)을 계산합니다.

  4. \(\mu \underline{\text{임금}}^{(g)} = \tilde{\beta_0}^{(g)} + \tilde{\beta_1}^{(g)}\underline{\text{임금}} + \tilde{\beta_2}^{(g)}\text{교육수준}^{(g)} + \tilde{\beta_3}^{(g)}(\underline{\text{임금}}\times\text{교육수준}^{(g)})\)을 계산합니다.

  5. \(\mu \overline{\text{임금}}^{(g)}\)과 \(\mu \underline{\text{임금}}^{(g)}\) 값을 저장해줍니다.

자, 그러면 이제 두 개의 벡터, \(\mu \overline{\text{임금}}^{(g)}\)과 \(\mu \underline{\text{임금}}^{(g)}\)를 가지게 되었습니다. 이 두 벡터를 요약해서 보여주기만 하면 됩니다. 그러면 우리는 임금이 매우 낮을 때와 임금이 매우 높을 때의 범죄율의 차이를 분포를 통해 직관적으로 확인할 수 있습니다. 이 방법이 효율적인 이유는 위의 4번과 5번에서처럼 \(\tilde{\beta}\)의 \(G\)개의 추출 결과를 가지고 그냥 \(\tilde{\textbf{b}}\)의 행렬과 결합한 뒤 계산만 해주면 되기 때문입니다. 아니면 그냥 바로 구한 \(\tilde{\beta}\) 행렬을 요약해서 보여줘도 되구요.

King et al. (2000)

이 파트에서는 King et al. (2000)의 주장과 제안을 간단하게 정리 및 소개하도록 하겠습니다. 나아가 그들이 (1) 왜 이런 논문을 썼는지, (2) 예측값(predicted values)와 기대값(expected values)의 차이점이 무엇인지, 그리고 1차 차분(first difference)이란 무엇인지, (3) 베이지안 접근법에 대해서 King et al. (2000)이 무어라 말하고 있는지 등에 대해 살펴보겠습니다.

library(mvtnorm)
model1 <- lm(log2(wdi_gdpcapcon2010) ~ wdi_trade + p_polity2 + 
               I(wdi_trade * p_polity2) + wdi_araland + wdi_agedr + 
               I(wdi_araland * wdi_agedr), data=QOG.s)
model1 %>% broom::tidy() %>% 
  mutate_if(is.numeric, round, 3)

set.seed(19891224)
beta_draws <- rmvnorm(n=1000, mean = coef(model1), sigma=vcov(model1))
head(beta_draws)

##      (Intercept)     wdi_trade    p_polity2 I(wdi_trade * p_polity2)
## [1,]    17.35005 -0.0023963830  0.025691386             7.312002e-04
## [2,]    18.68324 -0.0032424866  0.020011753             1.091112e-03
## [3,]    17.14560  0.0005784861 -0.008779808             8.476709e-04
## [4,]    18.22074  0.0042763462  0.017145542             7.995185e-04
## [5,]    16.62469  0.0070268007  0.086336501            -9.546159e-05
## [6,]    16.62564  0.0049103073  0.115611645            -3.416098e-05
##      wdi_araland   wdi_agedr I(wdi_araland * wdi_agedr)
## [1,] -0.01313041 -0.08070444              -0.0001032105
## [2,] -0.06878070 -0.09948163               0.0005771853
## [3,] -0.05279643 -0.08140526               0.0006251953
## [4,] -0.09277086 -0.10189863               0.0012669075
## [5,] -0.01525257 -0.07759279              -0.0001900096
## [6,] -0.04774735 -0.08325537               0.0004425008

trade <- quantile(QOG.s$wdi_trade, c(0.15, 0.85))
trade

##       15%       85% 
##  43.20651 121.37311

democracy <- c(seq(-10, 10, by=1))
democracy

##  [1] -10  -9  -8  -7  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6   7   8
## [20]   9  10

# 낮은 수준의 무역개방성
LowTrade <- cbind(1, trade[1], 
                  democracy, 
                  I(trade[1] * democracy),
                  mean(QOG.s$wdi_araland), mean(QOG.s$wdi_agedr), 
                  I(mean(QOG.s$wdi_araland) * mean(QOG.s$wdi_agedr)))
# 높은 수준의 무역개방성
HighTrade <- cbind(1, trade[2], 
                democracy, 
                I(trade[2] * democracy),
                mean(QOG.s$wdi_araland), mean(QOG.s$wdi_agedr), 
                I(mean(QOG.s$wdi_araland) * mean(QOG.s$wdi_agedr)))

LowTrade.ME <- t(LowTrade %*% t(beta_draws)) ## ME는 Marginal Effect입니다.
LT.mean <- apply(LowTrade.ME, 2, mean) ## 구해진 1,000개의 예측값의 평균
LT.se <- apply(LowTrade.ME,2, 
               quantile, c(0.025, 0.975)) ## 구해진 1,000개의 예측값의 표준편차

HighTrade.ME <- t(HighTrade %*% t(beta_draws))
HT.mean <- apply(HighTrade.ME, 2, mean)
HT.se <- apply(HighTrade.ME, 2, 
               quantile, c(0.025, 0.975)) ## 2.5/97.5 perentile

LT <- data.frame(Democracy=democracy,
                 Group = "Low Trade", ## 시뮬레이션이니 직접 하위 5%, 상위 5%의
                 # 관측치를 95% 신뢰구간을 위해 사용할 수 있습니다.
                 Mean=LT.mean, Lower=LT.se[1,], Upper= LT.se[2,])
HT <- data.frame(Democracy=democracy,
                 Group = "High Trade", ## 시뮬레이션이니 직접 하위 5%, 상위 5%의
                 # 관측치를 95% 신뢰구간을 위해 사용할 수 있습니다.
                 Mean=HT.mean, Lower=HT.se[1,], Upper= HT.se[2,])
Trade <- bind_rows(LT, HT)

Trade %>% 
  ggplot(aes(x=Democracy, y=Mean, color=Group, shape=Group)) + 
  geom_point() +
  geom_pointrange(aes(y = Mean, ymin = Lower, ymax = Upper)) + 
  scale_x_continuous(breaks = democracy)+
  scale_color_manual(values=futurevisions("mars")) +
  theme(axis.text.x  = element_text(vjust=0.5)) + 
  labs(x="민주주의 수준", y="경제 규모 (로그화)",
       caption="점추정치를 둘러싼 수직선은 95% 수준의 신뢰구간을 나타냄") +
  theme(legend.position = "bottom")

age <- quantile(QOG.s$wdi_agedr, c(0.15, 0.85))
age

##      15%      85% 
## 44.86155 82.48158

summary(QOG.s$wdi_araland)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1813  5.1117 12.0816 16.2298 22.6387 59.6467

agriland <- c(seq(0, 60, by = 5))


# 여기서 중요한 점은 예측변수의 행렬을 만들 때, 나중에 행렬곱셈을 해줄 시뮬레이
# 션된 계수값들의 순서는 OLS 분석에 투입된 변수 순서와 같다는 것입니다.
# 이를 고려해서 변수를 조작해주어야 합니다.

LowAge <- cbind(1, mean(QOG.s$wdi_trade), mean(QOG.s$p_polity2), 
                I(mean(QOG.s$wdi_trade) * mean(QOG.s$p_polity2)),
                agriland, age[1], 
                I(agriland * age[1]))

HighAge <- cbind(1, mean(QOG.s$wdi_trade), mean(QOG.s$p_polity2), 
                I(mean(QOG.s$wdi_trade) * mean(QOG.s$p_polity2)),
                agriland, age[2], 
                I(agriland * age[2]))

LowAge.ME <- t(LowAge %*% t(beta_draws))
LA.mean <- apply(LowAge.ME, 2, mean)
LA.se <- apply(LowAge.ME, 2, 
                   quantile, c(0.025, 0.975))

HighAge.ME <- t(HighAge %*% t(beta_draws))
HA.mean <- apply(HighAge.ME, 2, mean)
HA.se <- apply(HighAge.ME, 2, 
                   quantile, c(0.025, 0.975))

LA <- data.frame(Arable=agriland,
                 Group = "Low Age Dependency",
                 Mean=LA.mean, Lower=LA.se[1,], Upper= LA.se[2,])
HA <- data.frame(Arable=agriland,
                 Group = "High Age Dependency",
                 Mean=HA.mean, Lower=HA.se[1,], Upper= HA.se[2,])
Age <- bind_rows(LA, HA)

Age %>% 
  ggplot(aes(x=Arable, y=Mean, color=Group, shape=Group)) + 
  geom_point() +
  geom_pointrange(aes(y = Mean, ymin = Lower, ymax = Upper)) + 
  scale_x_continuous(breaks = agriland) +
  scale_color_manual(values=futurevisions("mars")) +
  theme(axis.text.x  = element_text(vjust=0.5)) + 
  labs(x="농지 (% 전체 토지)", y="경제 규모 (로그화)",
    caption="점추정치를 둘러싼 수직선은 95% 수준의 신뢰구간을 나타냄") +
  theme(legend.position = "bottom",
        legend.text=element_text(size=8.5))

Hainmueller et al. (2019)의 주장과 Brambor et al. (2006)에 대한 비판

Hainmueller et al. (2019)은 Brambor et al. (2006)가 곱셈을 통해 나타나는 상호작용항을 탐색하고 해석하는데 있어서 일종의 가이드라인을 제공하고는 있지만 몇 가지 중요한 이슈들을 간과하거나 언급조차 하고 있지 않다고 비판합니다. Hainmueller et al. (2019)에 따르면 그 문제는 크게 두 가지로 대별할 수 있습니다. 첫째, 선형 상호작용(linear interaction effect; LIE)에 대한 가정과 둘째, 충분한 정보량의 결여에 관한 것입니다.

Brambor et al. (2006)는 편미분을 취하는 방식을 통해서 상호작용항의 한계효과를 살펴보고 있습니다. 아래는 Hainmueller et al. (2019: 166)가 제시하고 있는 수식들로 고전적 선형회귀분석 모델에 곱셈 형태의 상호작용항이 포함된 모델을 보여줍니다.

\[\begin{equation} Y = \mu + \eta X + \alpha D + \beta(D\cdot X) + Z\gamma + \epsilon \end{equation}\]

이 모델에서 \(Y\)는 종속변수이고, \(D\)는 우리가 관심을 가지고 있는 핵심적인 예측변수, 혹은 처치변수입니다. \(X\)는 일종의 매개변수이고, \((D\cdot X)\)는 상호작용변수, \(Z\)는 일련의 통제변수들이라고 하겠습니다. \(\mu\), \(\epsilon\)은 각각 상수와 오차항을 보여줍니다. 이때, 핵심적인 예측변수 \(D\)의 종속변수 \(Y\)에 대한 한계효과는 다음과 같이 나타낼 수 있습니다.

\[\begin{equation} \text{ME}_D = \frac{\partial Y}{\partial D} = \alpha + \beta X \end{equation}\]

이 지점에서 Hainmueller et al. (2019)는 Brambor et al. (2006)이 간과한 점이 있다고 지적합니다. Brambor et al. (2006)의 논의에 따르면 우리는 단지 핵심 변수들의 한계효과가 일종의 선형 함수적 형태로 나타나리라고 가정해야 합니다. 그러나 Hainmueller et al. (2019)는 그 선형상호작용에 대한 가정은 결코 선험적인 것이 아니며 종종 유지되지도 않는다고 주장합니다. 따라서 Hainmueller et al. (2019)는 연구자들이 그들의 데이터를 한 번 더 살펴보고 한계효과를 진정 선형함수의 형태로 나타낼 수 있는지를 의심해보라고 주문합니다. Brambor et al. (2006)이 묵시적으로 LIE를 가정하고 많은 학자들이 그 가이드라인을 그저 따른다고 할지라도 Hainmueller et al. (2019)는 LIE는 결코 선험적으로 정당화될 수 없는 가정이며, 이를 확인하기 위해서는 연구자가 가진 데이터를 더 깊이 들여다보고 이해하는 것이 필요하다고 지적합니다.

이이서 Hainmueller et al. (2019)는 충분한 정보량의 문제에 대해서 지적합니다. 이 문제는 상호작용항의 효과를 살펴볼 수 있는 데이터의 범주 전반에 걸쳐서 실상 우리가 관측할 수 있는 가용성의 문제와 직결됩니다. 다르게 표현하자면, 만약 우리가 매우 치우친 형태의 분포를 가진 데이터를 가지고 있다면 그 데이터에서 매개변수의 한 값에서 우리는 한계효과가 명확하게 존재한다고 할만한 충분한 정보를 얻지 못할 수도 있기 때문입니다. 수학적, 행렬적 계산으로 도출해서 예측값의 변화를 보여줄 수는 있지만 과연 그것이 실제로는 존재하지 않는 데이터의 구간을 수리적 계산으로 그릴 뿐이라면? 한계효과에 대한 주장의 타당성에 의문을 제기할 수 있다는 것입니다. Hainmueller et al. (2019)는 논문 165쪽에서 두 가지 조건에 대해 명시하고 있습니다: “(1) 주어진 매개변수의 값에 대해서 \(X\) 값이 충분한 수의 관측치를 가지고 있어야 하며, (2) 그 매개변수의 값에서 핵심적인 예측변수, \(D\)의 변화가 존재해야 한다는 것”입니다. 특히 때로 매우 치우치거나 값이 일정하게 분포되어 있지 않은 자료를 사용하는 사회과학연구에서 이같은 노력을 주문하고 있습니다.

만약 주어진 매개변수의 특정한 값에 실질적으로 핵심적인 예측변수의 관측치들이 없거나, 거의 존재하지 않는다면 우리는 충분한 정보없이 한계효과를 그저 추정하는 것이 되고, 이를 Hainmueller et al. (2019: 165)는 “함수적 형태의 외삽 또는 내삽의 문제”라고 언급하고 있습니다. 정리하자면, Hainmueller et al. (2019)의 두 가지 핵심적인 주장은 첫째, Brambor et al. (2006)이 간과하거나 묵시적으로 가정하는 LIE의 문제를 고려할 것, (2) 상호작용 모델에 관한 이론적 이해와 우리가 사용하는 경험적 데이터 간의 간극을 좁힐 것으로 요약할 수 있습니다.

Hainmueller et al. (2019)의 대안