Introduction

주어진 x변수로 y값을 예측하는 경우를 생각해 보자. 이럴경우 y가 연속하는 값을 가지고 정규분포를 따른다면 가장 간단한단 모형으로 선형회귀모형을 생각할 수 있을 것이다. 그런데 그렇지 않고 y가 0, 1 두가지 값을 가지는 경우라면 어떻게 하겠는가? 생각할 수 있는 가장 간단한 방법은 $E(y) = P(y=1) = \beta_{0} + \beta_{1}x$ 이렇게 모형을 만드는 걸 생각할 수 있을것이다.

그런데 이렇게 모형을 만들 경우 생기는 문제점은 일단 $P(y=1)$ 이 값은 0에서 1사이의 값을 가지지만 $\beta_{0} + \beta_{1}x$ 이 값은 실수 전체값을 가질 수 있어서 극단적으로 확률값이 1보다 크거나 0보다 작은 모순적인 상황이 발생할 수 있다. 이러한 문제를 다루기 위해서 우리는 일반화 선형 모형(generalized linear model)을 알아볼 것이다.

일반화 선형 모형

구성 요소

본격적으로 들어가기에 앞서서 일반화 선형 모형을 구성하는 요소들에 대해서 알아보도록 하자.

  1. 확률 요소(random component)

독립적(independent)인 관측값(response variable) $(y_{1}, y_{2},…,y_{n})$ 들을 의미한다. 이때 이들의 분포는 특정 지수족을 따르고 동일한 분포를 따를 필요는 없다(non identical).

여기서 지수족을 따른 다는 말은 각 $y_{i}$ 가

\[\begin{aligned} f(y_{i}; \theta_{i}, \phi) = \exp \left(\frac{y_{i}\theta_{i} - b(\theta_{i})}{a(\phi)} +c(y_{i}, \phi) \right) \end{aligned}\]

이런 분포를 따른다는것을 의미한다.

여기서 $\phi$는 dispersion parameter라고 하며 이는 분산과 관련이 있다.

그리고 $\theta_{i}$는 natural parameter라고 하며 여기서 각 $y_{i}$마다 다름을 확인하자. 즉 $(y_{1}, y_{2},…,y_{n})$는 독립지만 동일한 분포를 따르진 않는다는것을 명심하자.

  1. 시스템 요소(systematic component) 혹은 선형 예측자(linear predictor)

예측변수(explanatory variable)과 그에 의한 선형함수를 의미한다. 선형회귀모형에서 $X\beta$를 의미한다 보면 된다. 이때 X는 확률변수가 아니라 상수라고 가정한다.

  1. 연결 함수(link function)

확률요소와 시스템 요소를 연결해주는 함수이다. 확률요소의 평균의 함수로써 이것을 시스템요소 같다고 둔다.

여기서 $g(u_{i}) = \theta_{i}$를 만족하는 연결 함수를 canonical link function이라하며 주로 이것을 사용해서 문제를 해결한다.

또한 연결 함수는 미분가능한 단조함수이다.

각 요소들의 예를 생각하면 선형 회귀 모형에서 우리가 일반적으로 $y_{i} = x_{i}’\beta + \epsilon_{i} \quad \epsilon_{i} \sim N(0, \sigma^2)$ 이렇게 모형을 만들었을때, 이것은 $E(y_{i}) = u_{i} =x_{i}’\beta$ 로 볼 수 있으므로 $g(u)=u$이고 시스템 요소는 $x_{i}’\beta$이다. 또한, $y_{i}$의 분포를 보면

\[\begin{aligned} f(y;u, \sigma) = \exp \left (\frac{uy - \frac{u^{2}}{2}}{\sigma^{2}} \right) - \frac{y^{2}}{2\sigma^{2}} - \frac{1}{2}log(2\pi\sigma^{2}) \end{aligned}\]

이다.

여기서,

\[\begin{aligned} &\theta = u \\ &b(\theta) = b(u) = u^{2} /2 \\ &\phi = \sigma^{2} \Rightarrow a(\sigma^{2}) = \sigma^{2} \end{aligned}\]

이 성립한다.

특히 중요히 봐야할 것은, $\theta\,,b(\theta)$ 이 둘을 유심히 보면 좋을것 같다.

확률 요소(random component)의 성질

다음으로 넘어가기 전에 확률 요소 $y_{i}$와 그 지수족 분포의 성질에 대해서 조금만 생각해보고 가도록하자. 여기서 각 지수족 분포의 점수 함수(score function)에 대해서 생각해보면 $L_{i} = \log f(y_{i};\theta,\phi) =\frac{y_{i}\theta_{i} - b(\theta_{i})}{a(\phi)} +c(y_{i}, \phi)$ 이므로,

\[\begin{aligned} & \frac{\partial \log f(y_{i};\theta,\phi)}{\partial \theta_{i}} = \frac{y_{i} - b'(\theta_{i})}{a(\phi)} \\ &\frac{\partial^{2} \log f(y_{i};\theta,\phi)}{\partial \theta_{i}^{2}} =\frac{- b''(\theta_{i})}{a(\phi)} \end{aligned}\]

이다. 또한 수리통계 내용에서 점수확률의 평균은 0이고 분산은 정보량이 됨이 알려져 있으므로, 다음이 성립한다

\[\begin{aligned} & E \left( \frac{\partial \log f(y_{i};\theta,\phi)}{\partial \theta_{i}} \right) = 0 \\ & E \left( \left\{ \frac{\partial \log f(y_{i};\theta,\phi)}{\partial \theta_{i}} \right\}^{2} \right) = E \left(-\frac{\partial^{2} \log f(y_{i};\theta,\phi)}{\partial \theta_{i}^{2}} \right) \end{aligned}\]

이를 통해서 다음이 성립함을 알 수 있다.

\[u_{i} = b'(\theta_{i})\] \[Var(Y_{i}) = a(\phi)b''(\theta_{i})\]

이것이 왜 필요하냐고 의문이 들 수 있는데 그것에 대해서 설명을 하자면 이로부터 해당 모델의 canonical link function을 유도가능하다. 위의 식의 첫번째 내용으로 부터 $b’^{-1}(u_{i}) = \theta_{i}$ 이므로 $g = b’^{-1}$가 이 모델의 cononical link function이 된다.

모델 설정

각 i번째 관측치에 대해서 설명변수 $x_{i} = (x_{i1}, x_{i2}, …,x_{ip})^{t}$ 가 주어졌다고 하자. 그리고 $\beta = (\beta_{1}, \beta_{2},…,\beta_{p})$가 추정해야 하는 모수이다.

이때, 일반화 선형 모형은

\[\begin{aligned} g(u_{i}) = \eta_{i} = x_{i}'\beta, \qquad i=1,...,n \end{aligned}\]

혹은 행렬 표현법에 의해서

\[g(u) = \eta = X\beta\]

이렇게 주어진다. 여기서 $X = (x_{1}, x_{2},..,x_{n})^T$는 $n * p$ 모델 행렬이다.

이 모형에서 주의해야 할점은 우변에 오차항이 붙지 않는다는 점이다. 좌변이 평균 $u$에 대한 함수로서 고정된 값이므로 우변 또한 고정된 값이 되어야하므로 확률변수인 오차항이 있으면 안된다.

최대가능도 추정(maximum likelyhood estimate)

이제 이렇게 만든 모델에서 우리는 $\beta$를 추정하는게 목적이다. 여기서 각 $y_{i}$들의 분포가 주어졌고 서로 독립이므로 그 결합확률분포를 구할 수 있으므로 최대가능도 추정법을 사용해서 $\beta$를 추정할것이다.

그러기 위해서 가능도함수를 최대화하면 되는데 여기서 지수족의 곱의 형태이므로 로그를 취한 로그가능도함수를 구해서 그것을 최대로 하는 $\beta$를 찾으면 된다. 로그가능도함수는 다음과 같다.

\[\begin{aligned} L(\beta) = \sum_{i=1}^{n} L_{i} = \sum_{i=1}^{n} \log f(y_{i};\theta,\phi) = \sum_{i=1}^{n} \left\{ \frac{y_{i}\theta_{i} - b(\theta_{i})}{a(\phi)} +c(y_{i}, \phi) \right\} \end{aligned}\]

여기서 이제 $L(\beta)$의 gradient가 0이되는 $beta$를 찾기 위해서는 각각의 $\beta_{j}$에 관해서, $\frac{\partial L}{\partial \beta_{j}}$를 구해야 하는데 $L(\beta)$의 식을 보면 알 수 있듯이 이식은 $\beta_{j}$에 대해서 바로 들어난 식이 아니다. 즉 저걸 구하기 위해서는 연쇄법칙(chain rule)을 사용해야한다.

이쯤 우리에게 주어진게 무엇이 있는지 확인하고 가보자.

\[\begin{aligned} & L(\beta) = \sum_{i=1}^{n} \left\{ \frac{y_{i}\theta_{i} - b(\theta_{i})}{a(\phi)} +c(y_{i}, \phi) \right\} \\ & u_{i} = b'(\theta_{i}) \\ & g(u_{i}) = \eta_{i} \\ & \eta_{i} = x_{i}^{t}\beta \end{aligned}\]

이 식들이 우리에게 주어진 조건들이다. 여기에 연쇄법칙을 적용하면,

\[\begin{aligned} \frac{\partial L}{\partial \beta_{j}} = \frac{\partial L}{\partial \theta_{i}} \frac{\partial \theta_{i}}{\partial u_{i}} \frac{\partial u_{i}}{\partial \eta_{i}} \frac{\partial \eta_{i}}{\partial \beta_{j}} \, \qquad j=1,2...,p \end{aligned}\]

이 성립한다. 여기서 index를 주의하도록하자. $i$는 1부터 n까지 가는 표본을 나타내는 index이고 $j$는 1부터 p까지 가는 특성의 개수를 나타내는 index이다. 이 둘을 잘 구분하도록 하자.

위의 것들을 정리하고 마지막으로 $Var(Y_{i}) = a(\phi)b’’(\theta_{i})$ 을 사용하면 다음과 같은 정규방정식(normal equation)을 구할 수 있다.

\[\begin{aligned} \sum_{i=1}^{n} \left\{ \frac{(y_{i} - u_{i})x_{ij}}{Var(Y_{i})} \frac{\partial u_{i}}{\partial \eta_{i}} \right\}=0 \, \qquad j=1,2...,p \end{aligned}\]

문제점은 이게 닫힌 형태 (closed form)이 존재하지 않아서 $\beta$를 명시적으로 구할 수가 없다는 없다는 점이다.

여기서 이 문제를 해결하는 세가지 방법이 알려져 있다.

  1. 뉴턴-랩손 방법(Newton-Raphson Algorithm)

  2. Fisher scoring method.

  3. 반복재가중최소제곱법.(Iteratiely reweighted least square, IRLS)

우리가 만약 모델을 만드는데 있어서 일반적인 연결함수가 아닌 canonical link function을 사용했을 경우에 세가지 방법이 모두 동치가됨이 알려져 있다. 다음절에 이를 증명해보자.

반복재가중최소제곱법.(Iteratiely reweighted least square, IRLS)

여기서 canonical link function을 사용했을 경우 뉴턴-랩손 방법이 Fisher scoring method가 되고 이를 풀었을때 반복재가중최소제곱법이 됨을 간략하게 설명할 것이다.

뉴턴-랩손 방법으로 $L(\beta)$가 최대가 되는 $\beta$를 얻기 위해서는

\[\begin{aligned} \hat{\beta}^{(m)} = \hat{\beta}^{(m-1)} - \left( \left. \frac{\partial^{2} L}{\partial \beta^{2}} \right|_{\beta=\hat{\beta}_{j-1}} \right)^{-1} \left(\left. \frac{\partial L}{\partial \beta} \right|_{\beta=\hat{\beta}_{j-1}} \right) \end{aligned}\]

을 반복적으로 시행하면 된다.

이 경우 편의상 \(U^{(m-1)} = \left( \left. \frac{\partial L}{\partial \beta} \right| _{\beta=\hat{\beta}_{j-1}} \right)\) 라고 하자.

또한 여기서 canonical link를 사용했을 경우 $\theta = \eta$가 되므로 $\frac{\partial L_{i}}{\partial \beta_{j}} = \frac{(y_{i} - u_{i})x_{ij}}{a(\phi)}$ 이 된다.

이로부터 $\frac{\partial^{2} L_{i}}{\partial \beta_{k} \partial \beta_{j}} = -\frac{x_{ik}x_{ij}}{a(\phi)} \frac{\partial u_{i}}{\partial \eta_{i}}$ 이 됨을 알 수있다.

이때, 이 값이 확률변수 $y_{i}$를 포함하지 않으므로 $E \left( \frac{\partial^{2} L_{i}}{\partial \beta_{k} \partial \beta_{j}} \right) = \frac{\partial^{2} L_{i}}{\partial \beta_{k} \partial \beta_{j}}$ 이 성립하므로 뉴턴-랩손 방법과 Fisher scoring method는 동치이다.

이제 Fisher scoring method와 반복재가중최소제곱법이 동치임을 유도해보자.

\[\begin{aligned} \hat{\beta}^{(m)} = \hat{\beta}^{(m-1)} + (I(\hat{\beta^{(m-1)}}))^{-1} \left(\left. \frac{\partial L}{\partial \beta} \right|_{\beta=\hat{\beta}_{j-1}} \right) \end{aligned}\]

를 반복해서 구해야한다.

이때 $I(\hat{\beta^{(m-1)}})$의 $(k, j)$ 번째 항을 보면

\[E\left(-\frac{\partial^{2} L_{i}}{\partial \beta_{k} \partial \beta_{j}}\right) = \frac{x_{ik}x_{ij}}{a(\phi)} \frac{\partial u_{i}}{\partial \eta_{i}} = \frac{x_{ik}x_{ij}}{Var(Y_{i})} \left(\frac{\partial u_{i}}{\partial \eta_{i}}\right)^2\]

이렇게 구해진다.

이 경우에, $w_{i} = \frac{1}{Var(Y_{i})} \left(\frac{\partial u_{i}}{\partial \eta_{i}}\right)^2$ 이고 $W = diag(w_{i})$라고 하면, $I(\beta) = X’WX$ 이렇게 쓸 수 있다.

\[\frac{\partial L_{i}}{\partial \beta_{j}} = \frac{(y_{i} - u_{i})x_{ij}}{Var(Y_{i})} \frac{\partial u_{i}}{\partial \eta_{i}}\]

이므로, $(y_{i} - u_{i}) \frac{\partial u_{i}}{\partial \eta_{i}} = v_{i}, \, v=(v_{1},v_{2},..,v_{n})$라 하면 $\frac{\partial L}{\partial \beta} = X’Wv$가 성립한다. 즉, 주어진 식이

\[\begin{aligned} &\hat{\beta}^{(m)} = \hat{\beta}^{(m-1)} + (X'WX)^{-1} X'Wv \\ &(X'WX)\hat{\beta}^{(m)} = (X'WX)\hat{\beta}^{(m-1)} + X'Wv \\ &X'WX\hat{\beta}^{(m)} = X'W(X\hat{\beta}^{(m-1)} + v) \\ &\hat{\beta}^{(m)} = (X'WX)^{-1}X'Wz \qquad (z=X\hat{\beta}^{(m-1)} + v) \end{aligned}\]

이렇게 되므로 Fisher scoring method와 반복재가중최소제곱법이 동치임을보였다.

예시

예시를 들기 위해서 다음과 같은 상황을 생각해보자.

어떤 한 도시에서 범죄 발생률을 연속된 3달간격으로 20번 측정했을때 첫번째 측정을 i라 하고 그때의 측정값이 1, 6, 16, 23, 27, 39, 31, 30, 43, 51, 63, 70, 88, 97, 91, 104, 110, 113, 149, 159 이렇게 주어졌다고 하자. 이때 $x_{i} = logi$로 하고 일반화 선형 모형을 적용시켜 보자.

여기서 반응변수가 범죄 발생률이므로 포아송 분포를 따른다고 예측할 수 있다. 평균이 $u$인 포아송 분포의 canonical link function은 $g(u)=\log u$임이 알려져 있다.

즉 우리는 $g(u_{i})=\log u_{i} = \beta_{0} + \beta_{1}x_{i}$라고 glm모델 만들것이다.

R 내장함수 사용

일단 먼저 R에 내장된 glm 함수를 통해서 $beta$를 추정하는 법을 알아보자.

y = c( 1, 6, 16, 23, 27, 39, 31, 30, 43, 51, 63, 70, 88, 97, 91, 104, 110, 113, 149, 159)
i = seq(1,20,1)
x = log(i)
glm_1 <- glm(y~x, family = poisson(link =log))
summary(glm_1)
## 
## Call:
## glm(formula = y ~ x, family = poisson(link = log))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0568  -0.8302  -0.3072   0.9279   1.7310  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.99600    0.16971   5.869 4.39e-09 ***
## x            1.32661    0.06463  20.525  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 677.264  on 19  degrees of freedom
## Residual deviance:  21.755  on 18  degrees of freedom
## AIC: 138.05
## 
## Number of Fisher Scoring iterations: 4

기본적으로 glm함수는 lm함수와 비슷하게 작동한다. 가장 큰 차이점은 family 매개변수인데 여기서 반응변수의 분포와 사용할 연결 함수를 지정할 수 있다. 우리는 포아송 회귀를 하고있기때문에 family = poisson(link =log)으로 진행하였다. 그 후 lm함수와 유사하게 summary함수를 사용해서 적합된 모형에 대한 정보를 알 수있다. 여기서 보면 절편항은 그 추정값이 0.996이고 $\beta_{1} = 1.32661$을 알 수 있다.

그 후 적합값 또한 구할 수있는데 이도 선형회귀와 마찬가지로 predict.glm함수를 사용해서 구할 수 있다. 이때 tpye매개 변수를 설정해 줄 수 있는데 type="response"는 $\hat{u}$을 계산하고 type="link"는 $g(\hat{u})$을 계산해준다. 기본값은 type="link"이다. 한번 확인해보자.

predict.glm(glm_1)                           # log(\hat{u})
##        1        2        3        4        5        6        7        8 
## 0.995998 1.915534 2.453428 2.835070 3.131094 3.372963 3.577461 3.754605 
##        9       10       11       12       13       14       15       16 
## 3.910857 4.050630 4.177069 4.292499 4.398685 4.496997 4.588524 4.674141 
##       17       18       19       20 
## 4.754566 4.830393 4.902119 4.970165
predict.glm(glm_1, type="link")             # log(\hat{u})$
##        1        2        3        4        5        6        7        8 
## 0.995998 1.915534 2.453428 2.835070 3.131094 3.372963 3.577461 3.754605 
##        9       10       11       12       13       14       15       16 
## 3.910857 4.050630 4.177069 4.292499 4.398685 4.496997 4.588524 4.674141 
##       17       18       19       20 
## 4.754566 4.830393 4.902119 4.970165
predict.glm(glm_1, type="response")          # \hat{u}
##          1          2          3          4          5          6          7 
##   2.707425   6.790563  11.628137  17.031585  22.899016  29.164829  35.782583 
##          8          9         10         11         12         13         14 
##  42.717356  49.941755  57.433612  65.174554  73.149058  81.343805  89.747218 
##         15         16         17         18         19         20 
##  98.349124 107.140501 116.113281 125.260201 134.574679 144.050720
log(predict.glm(glm_1, type="response"))    # log(\hat{u})
##        1        2        3        4        5        6        7        8 
## 0.995998 1.915534 2.453428 2.835070 3.131094 3.372963 3.577461 3.754605 
##        9       10       11       12       13       14       15       16 
## 3.910857 4.050630 4.177069 4.292499 4.398685 4.496997 4.588524 4.674141 
##       17       18       19       20 
## 4.754566 4.830393 4.902119 4.970165

R을 통해서 직접 구현하기

이번에는 이전에 말한 세가지 방법

  1. 뉴턴-랩손 방법(Newton-Raphson Algorithm)

  2. Fisher scoring method.

  3. 반복재가중최소제곱법.(Iteratiely reweighted least square, IRLS)

을 통해서 $\beta$를 직접구하고 실제로 그 값이 동일한지 확인해 볼것이다.

뉴턴-랩손 방법(Newton-Raphson Algorithm)

# 1. Newton-Raphson

X <- cbind(1, log(i))
beta_new <- c(1, 2) # 초기값 설
iter <- 1


while(iter<30){
    beta_old <- beta_new
    eta <- X %*% beta_old
    mu <- exp(eta)
    d.mu_d.eta <- exp(eta)
    var_y <- mu
    # w <- as.numeric((1/var_y) * (d.mu_d.eta)^2)
    W <- diag(as.vector(mu))
    z <- eta + (y-mu) * (1/d.mu_d.eta)
    beta_new <- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% z)
    crit = mean(abs(beta_old - beta_new))
    iter <- iter + 1
    if (crit <= 0.0001) {
        break
    }
}
beta_new
##          [,1]
## [1,] 0.995998
## [2,] 1.326610

Fisher scoring method.

#2 . Fisher scoring
X <- cbind(1, log(i))
beta_new <- c(1, 2)
iter = 1

while(iter < 30){
    beta_old <- beta_new
    eta <- X%*%beta_old
    mu <- exp(eta)
    d.mu_d.eta <- exp(eta)
    var_y <- mu
    w <- as.vector(mu)
    W <- diag(as.vector(w))
    v <- (y - mu) * (1/d.mu_d.eta)
    u <- t(X)%*%W%*%v
    I <- t(X)%*%W%*%X
    beta_new <- beta_old + (solve(I)%*%u)
    crit = mean(abs(beta_old - beta_new))
    iter <- iter + 1
    if (crit <= 0.0001) {
        break
    }
}
beta_new
##          [,1]
## [1,] 0.995998
## [2,] 1.326610

반복재가중최소제곱법.(Iteratiely reweighted least square, IRLS)

#3. IRLS

beta_new <- c(1,2)
iter = 1
while(iter < 30){
    beta_old <- beta_new
    eta <- X%*%beta_old
    mu <- exp(eta)
    d.mu_d.eta <- exp(eta)
    var_y <- mu
    w <- as.vector(mu)
    W <- diag(as.vector(w))
    z <- eta + (y - mu) * (1/d.mu_d.eta)
    beta_new <- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% z)
    crit = mean(abs(beta_old - beta_new))
    iter <- iter + 1
    if (crit <= 0.0001) {
        break
    }
}
beta_new
##          [,1]
## [1,] 0.995998
## [2,] 1.326610

실제로 위 세가지 방법으로 구한 $\hat{\beta}$값이 전부 일치하고, 실제로 glm함수를 통해서 구한 값과도 거의 유사함을 볼 수 있다.

Summary

오늘은 일반화 선형 모형을 사용하는 이유와 그것이 무엇인지 그리고 어떻게 추정하는지에 대해서 알아 보았다.

주로 많이 사용하게 되는건 분류문제(classification)에서 레이블(label)이 0과 1, 즉 이진으로 주어졌을때 로지스틱 회귀 (logistic regression)을 사용하게 될텐데 그건 단순히 glm에서 반응변수가 베루누이 분포를 따르고 연결함수를 $g(u) = \log (\frac{u}{1-u})$를 사용한 경우이다.

머신러닝 분야에서는 해당 문제를 풀때 확률적 경사 하강법(stochastic gradient descent)를 사용해서 푸는데 내 생각엔 그것 역시 결국에는 반복재가중최소제곱법을 조금 변형한거로 생각된다. 그 경우 학습률을 정하고 이전의 가중치(모수)에서 학습률*gradient를 뺀거로 해서 새롭게 가중치를 업데이트 하는데 이 역시 같은 결과를 이끌어 낸다.

이러한 이유에서 복잡한 딥러닝, 머신러닝 모델을 공부하기 전에 꼭 공부를 하고가면 도움이 될거라 생각한다.

댓글남기기