로지스틱회귀 MLE로 직접 풀어보기

지난 포스팅에서 최대가능도 추정법을 정규분포에 적용한 간단한 사례를 살펴보았습니다. 하지만 아무도 정규분포의 평균과 분산을 추정할 때 굳이 최대가능도 추정법을 사용하지 않을 것입니다. 표본평균, 표본분산이 있는데 왜 굳이 이런 힘든 수고를 겪어야 할까요? 게다가 이 값들은 최대가능도 추정치와 완전히 일치하거나, 조금만 조정하면 최대가능도 추정치로 만들 수 있는데 말입니다. 이런 경우에 최대가능도법을 굳이 적용하는 것은 닭 잡는 데 소 잡는 칼을 쓰는 격입니다.

사실 최대가능도 추정법은 이런 간단한 문제가 아닌, 표본에서 추정치를 직접 계산하기 어려운 좀 더 복잡한 사례에서 더 흔히 찾아볼 수 있습니다. 가장 흔히 접할 수 있는 사례는 일반화 선형모형 generalized linear model 입니다. 이 포스팅에서는 일반화 선형모형을 예로 들어 최대가능도 추정법을 어떻게 적용하는지 알아보겠습니다. 그리고 그 결과가 R의 내장함수가 산출하는 결과와 일치한다는 것까지도 확인해 보겠습니다.

일반화 선형모형에 대해 간단히 소개하자면, 이것은 일반적인 선형회귀분석의 가정이 충족되지 않는 경우에도 선형모형을 사용할 수 있게 해 줍니다. 이를테면 종속변수가 이분형 binary 이라든지, count data라든지 하는 경우에 말이죠. 로지스틱회귀 logistic regression 는 전자의 경우에 사용되는 가장 대표적인 일반화 선형모형입니다. 이 모형은 다음과 같이 표현할 수 있습니다:

p_i = P(Y_i = 1 | X_1 = x_1, …, X_p = x_p, beta) = 1 / (1 + exp(-(beta_0 + beta_1*x_1 + … + beta_p*x_p))) … (1)

여기서 Y_i는 이분형 종속변수 (0 또는 1의 값을 가짐), X_1부터 X_p까지는 p개의 예측(또는 독립)변수들을 나타냅니다. 로지스틱 회귀는 수식에서 알 수 있는 것과 같이 Y_i가 1의 값을 가질 확률을 독립변수들의 선형함수의 ‘역-로짓 변환’ inverse logit transformation 으로 예측하죠. 이 식을 좀 뒤틀어서 자료를 주어진 것으로, 모수치를 입력값으로 간주하면 가능도함수가 됩니다. 가능도함수는 따라서 다음과 같이 쓸 수 있습니다:

L(beta_0, …, beta_p | X, Y) = prod(p_i^Y_i*(1-p_i)^{1-Y_i}) … (2)

여기서 p_i는 (1)에서 계산한 바와 같습니다. (2)의 양변에 로그를 취하면 다음과 같이 변합니다:

log L(beta_0, …, beta_p | X, Y) = sum(Y_i*log(p_i)+(1-Y_i)*log(1-p_i)) … (3)

이것이 로그가능도 함수고, 회귀계수들에 대한 최대가능도 추정은 Y_i들이 주어졌을 때 우변을 최대화해 주는 beta들의 값을 찾음으로써 이루어집니다. 그런데 이 최대화 문제는 정규분포의 경우처럼 간단히 손으로 풀어서 계산하기 매우 힘들기 때문에, 컴퓨터를 동원할 수밖에 없습니다.

이제 R에서 이 문제를 풀어봅시다. 예제로 (또!) iris 데이터셋을 사용합니다. 반응변수를 이분형으로 만들기 위해 setosa 종의 자료는 제외하고, 나머지 100개의 데이터만 가지고 해 보겠습니다. 여기서 반응변수는 “iris의 종이 versicolor냐 아니냐”가 되겠습니다. 이것을 is_versicolor라는 변수에 저장합니다.

data <- iris[51:150,]
data$is_versicolor <- as.numeric(data$Species == ‘versicolor’)

이제 로그-가능도함수를 짜 봅시다. R의 optim 함수는 함수의 최소화만 가능하기 때문에, 로그가능도함수에 음의 부호를 취해야 한다는 것은 앞에서 언급했습니다:

neg_loglik <- function(data, parms){
Xb <- parms[1] + parms[2]*data[,1] + parms[3]*data[,2] + parms[4]*data[,3] + parms[5]*data[,4]
p <- 1 / (1+exp(-Xb))
-2*sum(dbinom(data$is_versicolor, 1, p, log=T))
}

로그-가능도함수의 정의에 맞게 함수를 작성했습니다. 마지막 줄에 -2를 곱해놓은 걸 보셨을 텐데, 여기에는 나중에 언급할 이유가 있습니다. (2는 곱하든 말든 최대화 문제에 영향을 주지 않습니다.) 이제 optim 함수를 사용해서 이 마이너스-로그-가능도 함수를 최소화 (로그-가능도함수를 최대화)해 봅시다:

init <- rep(0, 5)
fit <- optim(init, neg_loglik, data=data, control=list(maxit=10000))

이 최적화는 조금 오래 걸리기 때문에 최대 반복 횟수를 10,000까지 조금 늘려 주었습니다. 결과는 다음과 같습니다:

fit$par
[1] 42.599709 2.470021 6.675285 -9.429527 -18.270486

이 값들이 회귀계수들의 최대가능도 추정치입니다. 최소화된 -2 곱하기 로그가능도 값 (최대화된 로그가능도 값)은 다음과 같습니다:

fit$value
[1] 11.89856

이제 이 값들이 맞는지 R의 내장 함수인 glm을 사용해서 크로스체크해 봅시다. 결과는 다음과 같습니다:

fit_glm <- glm(is_versicolor ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data=data, family=’binomial’)
coef(fit_glm)
(Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
42.637804 2.465220 6.680887 -9.429385 -18.286137

완전하게는 아니지만 우리가 optim 함수를 써서 찾은 값들과 대체로 일치하는 것을 볼 수 있습니다. 다음 줄을 실행하면 deviance라는 값들이 마지막에 표시되는 것을 볼 수 있습니다:

summary(fit_glm)

해당 부분을 보면 다음과 같습니다:

Null deviance: 138.629 on 99 degrees of freedom
Residual deviance: 11.899 on 95 degrees of freedom

여기서 Residual deviance라는 값을 보면 아까 optim으로 최소화시킨 -2 곱하기 로그가능도 값과 같은 것을 확인할 수 있습니다. 이 값은 선형회귀에서 residual sum of squares와 같다고 보시면 됩니다. 즉 모델이 데이터를 제대로 예측하지 못한 정도를 나타냅니다. 구체적으로는 deviance는 앞에서 보신 것과 같이 -2 곱하기 로그가능도 값으로 정의됩니다.

참고로 Null deviance는 아무런 예측변수도 사용하지 않았을 때의 deviance 값입니다. (참고로 이런 통계 모형을 Null model이라고 합니다.) 아무런 예측변수도 없다면 가장 나이브한 예측 방식은 전 데이터셋의 비율을 그대로 활용하는 것입니다. 데이터셋에는 50개의 versicolor와 50개의 virginica가 있으므로, 나이브한 에측은 각각의 사례가 versicolor일 확률이 0.5라고 가정하는 것입니다. 이 값을 사용해서 deviance 값을 계산해 보면 다음과 같습니다:

-2*sum(dbinom(data$is_versicolor, 1, 0.5, log=T))
[1] 138.6294

위에서 본 Null deviance 값과 정확히 일치하는 것을 볼 수 있습니다. 그런데 deviance 값은 어떻게 해석할까요? 선형회귀의 경우와 같이, deviance 값의 감소는 예측되지 못한 부분의 감소, 즉 예측력의 개선으로 해석 가능합니다. 이 사례에서는 아무런 예측변수도 없을 때보다 4개의 예측변수를 썼을 때 deviance가 138.629 – 11.899= 126.73 만큼 개선되었다는 것을 의미합니다. (참고: Deviance 값이 로그가능도에 -2를 곱한 것으로 정의되는 이유는 이렇게 하면 deviance의 차이 값이 카이제곱 분포를 점근적으로 따른다는 것이 알려져 있기 때문입니다. 이 사실을 이용하여 가설검정을 할 수 있습니다.)

이와 같이 최대가능도법은 사실 어떤 마법 같은 것이 아니라 그냥 가능도함수를 최대화시키는 기법이라는 것을 알 수 있습니다. 그리고 자료와 가능도함수만 알면 최대화 과정은 컴퓨터에게 맡기고 누구나 구현 가능합니다. 참고로 로지스틱회귀의 경우에는 가능도함수가 최대화되는 회귀계수들의 값이 유일하게 결정될 수 있다는 것이 알려져 있습니다. 여기에는 몇 가지 조건이 있긴 하지만, 일상에서 맞닥뜨리는 문제들에는 대부분 해당되지 않기 때문에 안심하고 사용해도 됩니다.

About The Author

Related posts

Leave a Reply

Your email address will not be published. Required fields are marked *