R 에서 로지스틱 회귀분석 결과 보는 법
실습
내장데이터 turnout 데이터를 불러와보자.
turnout는 1992년 미국 총선에 대한 데이터로써, race(인종), 연령(age), 교육수준(educate), income(수입)에 따른 vote(투표여부)를 파악할 수 있다. 이 데이터는 투표를 했느냐 안 했느냐하는 종속변수에 관심이 있으므로 로지스틱 회귀분석을 사용할 수 있다.
로지스틱 회귀분석은 일반적인 회귀분석과 달리 glm() 함수를 통해 모형을 세운다. 이때 옵션으로써 family=binomial() 을 넣어주면 알아서 로지스틱 회귀분석을 해준다. 갑자기 ‘바이노미얼’이 나오는 이유는 종속변수가 0이냐 1이냐의 이항변수기 때문이다.
결과 해석
결과를 보면 일반적인 다중회귀분석과 비슷한듯 다른 점이 많다. 각 변수에 대한 회귀계수를 보는 법은 똑같지만 F검정이 없어졌고 이탈도deviance라는 것이 생겼다.
이탈도라는 것은 로지스틱 회귀모형이 얼마나 데이터를 못 설명하는지에 대한 척도로 보아도 무방하다. 영이탈도null deviance는 아무런 변수 없이 상수항만 있을 때의 이탈도로써 데이터가 전혀 없는 최악의 상황이라 할 수 있다. 이렇게 얻어진 잔차이탈도residual deviance는 작으면 작을수록 좋고, 카이제곱분포를 따르기 때문에 카이제곱 적합도 검정을 통해 모형이 적합한지 확인할 수 있다.
glm() 함수의 이름이 ‘일반화된 회귀분석’을 의미하긴 하지만, 기본적으로 회귀분석이기 때문에 다중공선성을 고려하지 않을 수 없다. 다행스럽게도 VIF는 vif() 함수를 사용함으로써 쉽게 구해낼 수 있다. 예제에선 다중공선성 문제가 없는 것으로 보인다.
모형 진단
반면 잔차그림이 매우 충격적인데, 다행스럽게도 로지스틱 회귀분석에서 잔차그림은 아무런 의미를 갖지 못한다. 따라서 굳이 확인해볼 필요도 없고, 확인해보았을 때 이상해도 신경쓸 필요가 없다. 앞서 언급한 적합도 검정이 있으면 충분하다.
이제 실제로 몇가지 방법을 통해 모형이 적합한지 확인해보자.
카이제곱 적합도 검정을 하는 방법은 아주 간단하다. 유의수준 $\alpha$ 에 대해 검정하고 싶다면 qchisq() 함수에 $(1 - \alpha)$ 와 잔차이탈도의 자유도를 넣어줘서 임계치를 계산하면 된다. 위와 같이 임계치보다 잔차이탈도가 작으면 모형은 적합하다고 본다.
호스머-렘쇼 적합도 검정은 로지스틱 회귀분석에서 사용하는 대표적인 적합도 검정으로, 마찬가지로 카이제곱 통계량을 통해 모델이 적합한지 검정해준다. ResourceSelection 패키지의 hoslem.test() 함수에 모형의 실제 종속변수와 적합치를 넣어줌으로써 검정할 수 있다. 귀무가설은 ‘모형이 적합하다’이므로 위와 같이 유의확률이 높아 귀무가설을 기각하지 못하면 모형이 적합한 것으로 본다. 그러나 프랭크 하렐에 따르면 여러가지 약점이 있어 이제는 더 이상 권장되지 않는다고 한다1.
rms 패키지의 lrm() 함수를 사용하면 모형의 우도비 검정을 해준다. 위 두가지 방법이 적합도 검정이고 귀무가설이 ‘모형이 적합하다’인 것과 달리 우도비 검정에선 귀무가설을 기각해야 모형이 적합한 것으로 본다. 이왕 rms 패키지를 쓰고 있다면 굳이 glm() 함수로 로지스틱 회귀모형을 만들고 그 아웃풋을 넣지 말고 그냥 lrm() 함수에 모형을 바로 넣어줘도 상관 없다. 예제로 보자면 lrm(vote~.,data=turnout) 을 입력함으로써 모든 과정을 한번에 수행해준다.
예측
마지막으로 로지스틱 회귀모형을 이용해 직접 확률을 계산해보도록 하자. 확률은 predict() 함수에 type='response' 옵션을 넣어줌으로써 간단하게 계산된다.
물론 우리는 이미 데이터를 갖고 있고 참값을 알고 있기 때문에 이러한 예측은 큰 의미가 없다. 따라서 newdata 옵션을 통해 새로운 데이터를 주고 그때도 잘 맞추는지 확인하는 과정이 필요하다.
코드
아래는 예제 코드다.
install.packages("Zelig")
install.packages("car")
install.packages("ResourceSelection")
install.packages('rms')
library(Zelig)
data(turnout); head(turnout)
out0<-glm(vote~.,family=binomial(),data=turnout); summary(out0)
library(car)
vif(out0)
win.graph(4,4); plot(out0$residuals,main="잔차그림")
qchisq(0.95,df=1995)
library(ResourceSelection)
hoslem.test(out0$y,fitted(out0))
library(rms)
lrm(out0)
lrm(vote~.,data=turnout)
predict(out0,type='response\')