로지스틱 회귀분석 조절효과 - lojiseutig hoegwibunseog jojeolhyogwa

[태그:] 로지스틱

종속변수가 0과 1로 된 이분변수일 때 Logit이나 Probit 같은 비선형모형 추정을 하게 된다. 그리고 연구의 관심이 특정 독립변수가 종속변수에 미치는 영향력이 또 다른 독립변수의 수준에 따라 어떻게 달라지는가 하는 것일 때 두 독립변수의 상호작용항을 설정하게 된다. 가장 대표적인 사례는 집단 변수와 시간 변수의 상호작용항이 설정되는 비선형 DID(difference-in-difference method, 이중차감법)이다. 그런데 비선형모형에서 상호작용 효과를 추정하는 것에 대해 많은 논란이 있다.

비선형모형에서 상호작용항에 대한 비판의 핵심은 “비선형모형에서 상호작용항 추정치는 방향, 크기, 통계적 유의도에서 실제 상호작용 효과를 평가하기에 부적절하다”는 것이며, 그 근거로는 Ai and Norton(2003)의 주장이 가장 널리 인용된다. Ai and Norton(2003)은 1980~1999년 사이에 JSTOR에 등재된 13개의 경제학술지에 게재된 논문 중 총 72개의 논문이 비선형모형에서 상호작용항을 설정하였는데, 적절한 해석을 한 경우는 하나도 없었다고 하였다. 이후에 발표된 DeLeire(2000)의 논문만이 예외라고 하였다.

Ai and Norton(2003)은 비선형모형에서 상호작용항 추정치가 부적절한 이유를 다음과 같은 네 가지로 제시하였다. 첫째, 상호작용항의 회귀계수가 0일 경우에도 상호작용 효과는 0이 아닐 수 있다. 둘째, 상호작용 효과의 통계적 유의도는 상호작용항 회귀계수에 대한 단순한 t-test로 검증할 수 없다. 셋째, 상호작용 효과는 선형모형에서와 달리 독립변수에 조건적(conditional)이다. 이는 비선형모형에서 상호작용항이 아닌 일반적인 독립변수의 한계효과(marginal effect)가 독립변수에 조건적인 것과 동일하다. 넷째, 상호작용 효과는 다른 공변인의 값에 따라 다른 부호를 가질 수 있다. 따라서 상호작용항 회귀계수의 방향이 반드시 상호작용 효과의 방향을 알려주는 것은 아니다.

Ai and Norton(2003)에 따르면 비선형모형에서 상호작용 효과를 구할 때 가장 흔한 실수는 상호작용항을 일반적인 변수와 동일하게 취급하여 미분으로 한계효과를 구하는 것이다. 이들은 실제 상호작용 효과는 교차미분(cross derivative 혹은 cross difference)으로 구해야 하며, 그 표준오차는 Delta method를 적용하여 추정할 것을 제안하였다. 그리고 이들은 이러한 방법에 의해 Stata에서 상호작용 효과를 구하는 inteff라는 명령어를 제안하였다(Norton et al., 2004).

Google Scholar에서 확인해보니 Ai and Norton(2003)의 논문은 총 3,122회, Web of Science에서만 1,116회가 인용되었다(2015.5.14. 기준). 실로 엄청난 영향력이라 할 만하다. 하지만 이후 많은 학자들이 이들의 주장을 보완 혹은 비판하거나 대안적인 방법을 제안하기도 하였다. 이 중 몇 가지를 살펴보자.

우선 DeLeire(2004)는 비선형 DID에서 한계효과로 상호작용 효과를 구할 때 Ai and Norton(2003)이 제안한 방법을 대부분 옹호하였다. 하지만 LPM(linear probability model, 선형확률모형)도 DID의 불편 추정치를 구하는 좋은 대안으로 제시하였다. 종속변수가 이항변수이면서 DID를 할 경우 많은 사람들이 LPM을 선택하곤 하는데, 이 때 DeLeire(2004)의 논문이 그 근거로 인용되기도 한다.

Ai and Norton(2003)의 주장을 가장 정면으로 반박한 논문으로는 Puhani(2008)가 있다. Puhani(2008)는 비선형 DID에서 Ai and Norton(2003)이 제안한 cross-difference는 실제 처치효과(treatment effects on the treated)가 아니며, 상호작용항 회귀계수의 방향은 처치효과의 방향과 일치한다고 주장하였다. 따라서 Logit, Probit, Tobit과 같은 비선형 DID에서 상호작용항의 회귀계수에 초점을 두는 것이 타당하다고 하였다. 다만 통계적 유의도는 별개의 문제로 Ai and Norton(2003)이 제안한 것과 같이 Delta method로 구하는 것이 좋다고 하였다. 이후 Karaca-Mandic et al.(2012)가 Puhani(2008)의 주장을 인정하였다. 이 논문의 공동 저자에는 Norton이 있었다. 이들은 Ai and Norton(2003)의 주장을 기본적으로 채택하였지만, 비선형 DID 모델은 특별한 경우로 이 때 처치효과는 상호작용항의 회귀계수가 대표한다는 것을 확인하였다. 다만 비선형 DID를 Logit으로 추정할 경우 상호작용을 odds ratio로 해석하는 것은 매우 복잡하기 때문에 추천하지 않았다.

DeLeire(2004), Puhani(2008), Karaca-Mandic et al.(2012)의 글은 대부분 비선형 DID에서 처치효과를 구하는 것에 초점을 두었다. 즉 확률차이로 해석할 수 있는 한계효과를 구하는 것의 문제인 것이다. 이에 대해서는 Ai and Norton(2003)의 주장에 대한 보완이나 반박이 주를 이루었는데, Blundell et al.(2004), Athey and Imbens(2006)는 비선형 DID에서 처치효과를 구하는 대안적인 방법을 제시하기도 하였다. (개인적으로 완전히 이해하지 못해 내용을 소개하지는 않는다.)

하지만 비선형 DID에서 처치효과라는 특별한 경우가 아니라 일반적으로 비선형모형에서 상호작용 효과를 어떻게 파악해야 하는가에 대해서는 Buis(2010)와 Greene(2010)의 글이 더욱 포괄적이고 풍부한 시사점을 제공하는 것 같다.

Ai and Norton(2003)은 비선형모형에서 상호작용 효과를 한계효과(marginal effect)로 해석하는 것에 초점을 두었지만, Buis(2010)는 비선형모형에서 상호작용효과는 additive한 marginal effect로 해석할 수도 있고 multiplicative한 odds ratio(승산비)로 해석할 수도 있다는 것을 강조하였다. 적당한 우리말이 생각나지 않아 그냥 영문으로 옮겼는데, 한계효과를 구하여 확률의 ‘차이’로 해석할 수도 있고 odds ratio를 구하여 승산의 ‘변화율(혹은 비)’로 해석할 수도 있다는 말이다. Ai and Norton(2003)이 언급한 바와 같이 상호작용항에 대한 한계효과와 odds ratio는 그 방향이 달라질 수 있는데, Buis(2010)는 두 가지 모두가 실제 상호작용 효과에 대한 정확한 대푯값이며 무엇을 보고할 것인지는 실질적인 연구질문에 달려있는 것이라 정리하였다.

Buis(2010)도 논문에서 사례를 들어 설명했는데, 그것보다는 더 간단한 예를 들어보자. 만약 종속변수가 고용 여부, 독립변수가 장애 여부, 그리고 조절변수로 취업서비스를 받았는지 여부가 있다고 하자. 그리고 취업서비스를 받지 않은 경우, 받은 경우에 비장애인과 장애인의 고용확률 p가 아래와 같다고 하자. 상호작용 효과를 확률차이로 계산하면 (0.5-0.4)-(0.9-0.8)=0.0이다. 즉, 취업서비스가 비장애인-장애인의 고용확률 격차에 아무런 영향을 주지 않는다.

비장애인 : 취업서비스 X p=0.8, 취업서비스 O p=0.9
장애인   : 취업서비스 X p=0.4, 취업서비스 O p=0.5

하지만 동일한 경우를 odds ratio로 해석해보자. 일단 odds는 아래와 같이 구할 수 있다. 그리고 취업서비스의 odds ratio는 비장애인의 경우 9/4=2.25, 장애인은 1/0.67=1.5이다. 즉, 비장애인은 취업서비스를 받았을 경우 받지 않은 경우에 비해 고용될 odds가 2.25배이고, 장애인은 1.5배이다. 그리고 상호작용 효과는 odds ratio의 ratio가 된다. 즉, 1.5/2.25=0.67이다. 따라서 취업서비스는 비장애인과 장애인의 상대적인 고용 승산 변화율에 부적인 영향을 준다.

비장애인 : 취업서비스 X odds=0.8/(1-0.8)=4, 취업서비스 O odds=0.9/(1-0.9)=9
장애인   : 취업서비스 X odds=0.4/(1-0.4)=0.67, 취업서비스 O odds=0.5/(1-0.5)=1

동일한 자료에서 상호작용 효과를 한계효과로 판단했을 때와 odds ratio로 판단했을 때 결과가 전혀 달라지는 것이다. 그리고 중요한 것은 Logit에서 상호작용항의 odds ratio가 갖는 의미는 ‘ratio of odds ratio’라는 것이다. 따라서 이를 ‘확률차이’로 해석해서는 안된다.

Greene(2010)은 이에 대해 더욱 포괄적인 관점을 제공한다. Greene(2010)은 Ai and Norton(2003)이 제안한 cross difference는 계산방법이 틀린 것은 아니지만 변수의 관계 측면에서 해석이 어렵다는 점을 지적하였다. 즉, ‘단위 변화’라는 partial effect의 의미를 부여할 수 없다는 것이다. Greene(2010)은 응용계량분석은 다음의 두 가지 단계로 진행할 것을 제안하였다.

1. 적절한 통계적 절차와 원칙에 의해 모델을 설정한다. 모델 특정화에 대한 통계적 검증은 이 단계에서 이루어지며, 가설 검증은 모델의 회귀계수와 모델 특정화의 구조적 측면에 대해 이루어진다. partial effect는 회귀계수도 아니고 모델 특정화의 요소도 아니다. 그것은 특정화되고 추정된 모델의 함의(implication)이다.

2. 모델을 추정하였다면 회귀계수값, 예측값, partial effect, 상호작용과 같은 모델의 함의를 독자에게 알린다. 이러한 목적을 위해서 통계적 수치에 부가되는 시각적 표현은 매우 유용하다. 가설검증은 이 지점에서 필요하지 않다. partial effect가 궁극적인 추정의 목적이라 할지라도, 모델 설계자가 모델에 의해 산출될 partial effect에 대한 가설을 설정하여 구조적 모델을 설계하는 경우는 거의 없다.

Greene(2010)은 본인의 이러한 주장이 대부분의 통계 소프트웨어가 추정된 partial effect에 대한 표준오차와 t값을 제시하고, 그것을 보고하는 일반적인 관행과 배치된다고 인정하였다. 하지만 이렇게 손쉽게 계산이 가능함에도 불구하고 분석에서 가설검정이 유용한 지점은 분석 단계가 아니라 모델 구성 단계에 있다고 하였다.

워낙 여러 주장이 대립되고 아직까지 완전한 합의가 없어 더욱 혼란스러울 수도 있겠다. 하지만 Ai and Norton(2003)을 비롯하여 한계효과에만 초점을 두었던 초기 주장과 달리, 한계효과와 odds ratio가 갖는 의미의 차이에 주목한 Buis(2010), 가설검증은 회귀계수에 대해 이루어져야 하고 한계효과나 odds ratio는 독자에게 어떠한 함의를 전달할 것인가에 대한 선택이라는 Greene(2010)의 주장이 더 설득력이 있어 보인다.

참고문헌

  • Ai, C. and Norton, E. C. (2003). Interaction term in logit and probit models. Economics Letters, 80(1), 123-129.
  • Athey, S. and Imbens, G. W. (2006). Identification and inference in nonlinear difference-in-differences models. Econometrica, 74(2), 431-497.
  • Buis, M. L. (2010). Stata tip 87: Interpretation of interactions in nonlinear models. The Stata Journal, 10(2), 305-308.
  • DeLeire, T. (2000). The wage and employment effects of the Americans with Disabilities Act. Journal of Human Resources, 35(4), 693-715.
  • DeLeire, T. (2004). A note on calculating difference in differences using probit models versus linear probability models. Michigan State University.
  • Greene, W. (2010). Testing hypotheses about interaction terms in nonlinear models. Economics Letters, 107(2), 291-296.
  • Karaca-Mandic, P., Norton, E. C. and Dowd, B. (2012). Interaction terms in nonlinear models. Health Services Research, 47(1pt1), 255-274.
  • Norton, E. C., Wang, H. and Ai, C. (2004). Computing interaction effects and standard errors in logit and probit models. The Stata Journal, 4(2), 154-167.
  • Puhani, P. A. (2008). The treatment effect, the cross difference, and the interaction term in nonlinear “difference-in-differences” models. IZA Discussion Paper No. 3478.

로지스틱 회귀분석을 배울 때 최적의 R2에 대한 합의가 없고 제시하지 않는 경우도 많다고 한 게 기억난다. 그래서 나도 로지스틱 회귀분석을 여러 차례 사용했지만 R2는 한번도 제시한 적이 없었다. 하지만 R2를 제시한 논문들도 많은데, 그 종류가 제각각이다. 가장 많이 볼 수 있는 것은 pseudo R2, Cox & Snell R2, Nagelkerke R2 세 가지인 것 같다. 이것들이 어떤 특성을 갖고 차이는 무엇인지 참고가 될 듯 하여 Paul D. Allison 교수(University of Pennsylvania)의 글을 소개한다.

http://www.statisticalhorizons.com/r2logistic

What’s the Best R-Squared for Logistic Regression?

February 13, 2013 By Paul Allison

로지스틱 회귀분석에 대해 내가 가장 많이 받은 질문 중의 하나는 “내 모형이 데이터에 적합(fit)한지 어떻게 말할 수 있는가?”이다. 이 질문에 대한 답변으로 두 가지의 일반적인 접근방식이 있다. 하나는 독립변수에 기반하여 종속변수를 얼마나 잘 예측하는지를 측정하는 것이다. 다른 하나는 모형이 더 복잡해질 필요가 있는지, 구체적으로 데이터를 잘 대표하기 위하여 추가적인 비선형항(nonlinearity)과 상호작용항이 필요한지를 검증하는 것이다.

이후 글에서 모형 적합에 대한 두번째 접근방식을 논의할 것이며, 왜 내가 Hosmer-Lemeshow 적합도 검증(goodness-of-fit test)을 좋아하지 않는지 설명할 것이다. 이 글에서는 예측력에 대한 R2 측정치에 초점을 둘 것이다. 그 과정에서 나는 이러한 측정치와 관련하여 오랫동안 추천해온 것을 철회할 것이다.

불행히도, 로지스틱 회귀분석에 대한 R2를 계산하는 많은 방법이 있는데 어떠한 것이 최선인지에 대한 합의는 없다. Mittlbock and Schemper(1996)는 12개의 측정치를 검토하였고, Menard(2000)도 몇 개를 검토하였다. 통계 소프트웨어에서 가장 많이 보고되는 방법은 McFadden(1974)이 제안한 것, 그리고 Cox and Snell(1989)이 제안한 것과 그것의 “수정된” 버전으로 보인다. 하지만 Cox-Snell R2(수정 버전과 비수정 버전 모두)는 실제로 그 이전에 Maddala(1983)와 Cragg and Uhler(1970)가 제안한 것이었다.

통계 소프트웨어 중에서 SAS와 Statistica는 Cox-Snell 측정치를 보고한다. JMP와 SYSTAT은 McFadden과 Cox-Snell 측정치를 모두 보고한다. SPSS는 이항 로지스틱 회귀분석에는 Cox-Snell 측정치를 보고하지만 다항로짓과 순위로짓에는 McFadden 측정치를 보고한다.

수년동안 나는 McFadden R2보다는 Cox and Snell R2를 추천해왔지만, 최근에는 그것이 실수였다고 결론내렸다. 나는 이제 McFadden R2가 더 좋은 선택이라고 믿는다. 하지만 나는 좋은 특성을 가지고, 직관적으로 와닿고, 쉽게 계산할 수 있는 다른 R2도 알게 되었다. 현 시점에서 나는 그것을 McFadden R2보다 좋아한다. 하지만 나는 그것에 대해 더 알게 되기 전까지는 확실한 추천을 하지는 않으려고 한다.

이제 구체적으로 살펴보자. 로지스틱 회귀분석은 물론 우도함수(likelihood function)를 최대화함으로써 추정된다. L0를 아무런 예측변수로 포함하지 않은 모형에 대한 우도함수의 값이라고 하고, LM을 추정한 모형에 대한 우도라고 하자. McFadden R2는 다음과 같이 정의된다.

R2McF = 1 – ln(LM) / ln(L0)

이 공식의 원리는 ln(L0)이 선형 회귀분석에서의 잔차제곱합(residual sum of squares)과 동일한 역할을 한다는 것이다. 결과적으로 이 공식은 “잔차 분산(error variance)”의 감소비율에 해당한다. 이 측정치는 때로는“유사(pseudo)” R2라고 불리기도 한다.

Cox and Snell R2는 다음과 같이 정의된다(n은 표본수를 의미한다).

R2C&S = 1 – (L0 / LM)2/n

이 공식의 원리는 선형 회귀분석과 동일하다. 다른 말로 하면, 선형 회귀분석에서의 일반적인 R2는 정확히 이 공식에 의해 예측변수가 없는 모형과 있는 모형의 우도로 결정된다는 것이다. 그렇다면 이 측정치는 유사(pseudo) R2보다는 “일반화(generalized)” R2라고 부르는 것이 적절할 것이다. 이와 대조적으로, McFadden R2는 OLS R2를 특별한 사례로 가지고 있지 않다. 나는 이러한 Cox-Snell R2의 특성이 매우 매력적이라고 생각했는데, 특히 이 공식이 최우(maximum likelihood) 추정을 하는 다른 종류의 회귀분석에 자연스럽게 확장될 수 있기 때문이다(예를 들어 가산자료(count data)에 대한 negative binomial regression이나 생존자료(survival data)에 대한 Weibull regression).

하지만 잘 알려져 있듯이 Cox-Snell R2는 상한값이 1.0보다 작다는 큰 문제가 있다. 상한값은 1 – L02/n이 된다. 이는 1.0보다 상당히 작을 수 있으며, 오직 사건이 발생한 사례의 비중인 p에 의존한다.

upper bound = 1 – [pp(1-p)(1-p)]2

p=0.5일 때 0.75의 최대값을 가지며, p=0.9 혹은 p=0.1일 경우 상한값은 0.48에 불과하다.

선형 모형에서의 R2와 같은 R2를 원했던 사람들에게 이것은 매우 실망스러운 것이다. 이를 단순하게 교정하는 방법이 있는데, Nagelkerke(1991)가 제안한 것으로 R2C&S를 그것의 상한값으로 나눠주는 것이다. 하지만 이러한 교정은 임시변통일 뿐이며 애초 R2C&S가 가진 이론적 장점을 심각하게 훼손한다. 또한 이는 일반적으로 선형확률모형(linear probability model)에서 얻을 수 있는 것과 비교하여 오해를 사기 쉬울 정도로 높다(하지만 어떤 이는 이를 하나의 특징으로 보기도 한다).

따라서 약간 주저되긴 했지만, 나는 McFadden 편으로 넘어가기로 결정했다. Menard(2000)가 주장했듯이, 그것은 Kvalseth(1985)가 제시한 좋은 R2의 8가지 기준을 대부분 충족한다. p=0.5일 때 McFadden R2는 수정하지 않은 Cox-Snell R2보다 약간 작은 경향이 있다. p가 0이나 1에 근접할 때 McFadden R2는 커지는 경향이 있다.

하지만 최근 Tjur(2009)가 제안한 새로운 R2가 있는데, 나는 이것을 McFadden R2보다 선호한다. 그것은 직관적인 장점이 많고, 상한값이 1.0이며, 선형 모형에서 R2 정의와 밀접히 관련된다. 또한 계산하기도 쉽다.

정의는 매우 단순하다. 종속변수의 두 가지 범주 각각에 대해 사건의 예측확률의 평균값을 계산한다. 그 다음에 두 평균값의 차이를 구한다. 그것이 전부다!

이것의 취지는 명료하다. 어떤 모형이 좋은 예측을 한다면, 사건이 발생한 사례들은 높은 예측값을 가져야 하고 사건이 발생하지 않은 사례들은 낮은 예측값을 가져야 한다. 또한 Tjur는 이 R2(그는 식별계수(coefficient of discrimination)라 불렀다)가 잔차제곱에 기반한 두 가지 R2 공식의 산술평균과 동일하며, 잔차제곱에 기반한 다른 두 가지 R2 공식의 기하평균과 동일하다는 것을 증명하였다.

다음은 Tjur의 R2를 Stata에서 계산한 사례이다. 나는 753명의 기혼여성의 경제활동참여에 대한 잘 알려진 데이터세트(Mroz, 1987)를 사용하였다. 종속변수인 inlf는 여성이 경제활동참여 상태이면 1, 아니면 0으로 코딩되어 있다. 로지스틱 회귀분석에는 여섯 개의 예측변수가 포함되었다. 명령문은 다음과 같다.

use "http://www.uam.es/personal_pdi/economicas/rsmanga/docs/mroz.dta", clear
logistic inlf kidslt6 age educ huswage city exper
predict yhat if e(sample)
ttest yhat, by(inlf)

predict 명령어는 예측값을 산출하여 그것을 새로운 변수인 yhat에 저장한다. if e(sample)은 회귀모형에 포함되지 않은 사례의 예측값을 계산하는 것을 방지한다. ttest 명령어는 두 집단의 예측값 평균의 차이를 구하는 가장 쉬운 방법이다(p-value는 무시하라). 경제활동 참여자의 예측값 평균은 0.680이고, 미참여자의 예측값 평균은 0.422로 나타난다. 그 차이인 0.258이 Tjur의 R2이다. 다른 R2와 비교하면, Cox-Snell R2는 0.248, McFadden R2는 0.208, 수정된 Cox-Snell R2는 0.332이다.

아래는 동일한 내용의 SAS 코드이다.

proc logistic data=my.mroz;
    model inlf(desc) = kidslt6 age educ huswage city exper;
    output out=a pred=yhat;
proc ttest data=a;
    class inlf; var yhat; run;

Tjur R2에 대한 가능한 반대 이유의 하나는 Cox-Snell R2, McFadden R2와 달리 최우함수의 수치에 기반하지 않는다는 것이다. 이로 인해 모형에 변수를 추가해도 Tjur R2는 줄어드는 것이 가능하다. 하지만 Kvalseth(1985)는 R2는 특정 추정방법에 기반하지 않는 것이 바람직하다고 하였다. 그렇다면 서로 다른 방법을 사용하여 예측을 하는 것이 모형의 예측력을 비교하는데 정당할 것이다. 예를 들어 로지스틱 회귀분석에 기반한 예측과 분류트리기법(classification tree method)에 의한 예측을 비교할 수 있다.

또다른 잠재적인 불만은 Tjur R2를 순위로짓이나 다항로짓에 쉽게 일반화할 수 없다는 것이다. McFadden R2와 Cox-Snell R2는 그러한 일반화가 간단하다.