logistic-regression.Rmd이분 로지스틱 회귀모형은 종속변수가 2가지 범주를 취하는 경우에 사용된다.
\(N\)개의 객체로 이루어진 학습데이터 \(\{(\mathbf{x}_i, y_i)\}_{i = 1, \cdots, N}\)를 아래와 같이 정의하자.
\(\mathbf{x}_i\) 관측값을 이용하여 \(y_i\)의 기대값 \(P_i\)을 추정하는 모형을 아래와 같이 로지스틱 함수로 정의하자.
\[\begin{eqnarray} P_i &=& P(y_i = 1 \,|\, \mathbf{x}_i)\\ &=& E[y_i | \mathbf{x}_i]\\ &=& \frac{\exp(\beta_0 + \boldsymbol\beta^\top \mathbf{x}_i)}{1 + \exp(\beta_0 + \boldsymbol\beta^\top \mathbf{x}_i)} \tag{1.1} \end{eqnarray}\]
여기에서 \(\boldsymbol\beta \in \mathbb{R}^{p}\)는 \(\mathbf{x}_i\)와 동일한 차원의 벡터이다 (\(\boldsymbol\beta = [\beta_1 \, \beta_2 \, \cdots \, \beta_p]^\top\)).
식 (1.1)는 모든 \(\mathbf{x}_i\)값에 대해 0에서 1 사이의 값을 갖게 되므로 각 범주에 속할 확률을 추정하는 데 적합한 반면, 변수 \(\mathbf{x}\) 및 계수들에 대해 선형이 아니므로 추정이 어렵다. 그러나 아래와 같이 로짓(logit) 변환을 통해 선형회귀식으로 변환할 수 있다.
\[\begin{eqnarray} logit(P_i) &=& \ln \left[ \frac{P_i}{1 - P_i} \right]\\ &=& \ln(\exp(\beta_0 + \boldsymbol\beta^\top \mathbf{x}_i))\\ &=& \beta_0 + \boldsymbol\beta^\top \mathbf{x}_i \tag{1.2} \end{eqnarray}\]
식 (1.2)에서 확률 \(P_i\)는 직접적으로 관측되는 것이 아니고 0 또는 1을 갖는 \(y_i\)가 관측되므로, \(P_i\)를 일종의 잠재변수(latent variable)로 해석할 수 있다.
\[\begin{equation} y_i = \begin{cases} 1 & \text{ if } \beta_0 + \boldsymbol\beta^\top \mathbf{x}_i + \varepsilon_i > 0 \\ 0 & \text{ otherwise } \end{cases} \tag{1.3} \end{equation}\]
식 (1.3)에서 \(\varepsilon_i\)는 표준로지스틱분포(standard logistic distribution)을 따른다.
data(student, package = "dmtr")| id | x1 | x2 | x3 | y |
|---|---|---|---|---|
| 1 | 0 | 8 | 2 | 우수 |
| 2 | 1 | 7 | 1 | 우수 |
| 3 | 0 | 9 | 0 | 우수 |
| 4 | 1 | 6 | 4 | 우수 |
| 5 | 1 | 8 | 2 | 우수 |
| 6 | 0 | 7 | 3 | 우수 |
| 7 | 0 | 7 | 0 | 보통 |
| 8 | 1 | 6 | 1 | 보통 |
| 9 | 0 | 7 | 2 | 보통 |
| 10 | 0 | 8 | 1 | 보통 |
| 11 | 0 | 5 | 2 | 보통 |
| 12 | 1 | 8 | 0 | 보통 |
| 13 | 0 | 6 | 3 | 보통 |
| 14 | 1 | 7 | 2 | 보통 |
| 15 | 0 | 6 | 1 | 보통 |
로지스틱 모형에서 회귀계수의 추정을 위해서 주로 최우추정법(maximum likelihood estimation)이 사용된다. \(N\)개의 객체로 이루어진 학습데이터에 대해 우도함수는 다음과 같다.
\[\begin{equation*} L = \prod_{i = 1}^{N} P_i^{y_i} (1 - P_i)^{1 - y_i} \end{equation*}\]
그리고 우도함수에 자연로그를 취하면 아래와 같이 전개된다.
\[\begin{eqnarray} \log L &=& \sum_{i = 1}^{N} y_i \log P_i + \sum_{i = 1}^{N} (1 - y_i) \log (1 - P_i)\\ &=& \sum_{i = 1}^{N} y_i \log \frac{P_i}{1 - P_i} + \sum_{i = 1}^{N} \log (1 - P_i)\\ &=& \sum_{i = 1}^{N} y_i (\beta_0 + \boldsymbol\beta^\top \mathbf{x}_i) - \sum_{i = 1}^{N} \log (1 + \exp (\beta_0 + \boldsymbol\beta^\top \mathbf{x}_i) )\\ &=& \sum_{i = 1}^{N} y_i \left(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij} \right) - \sum_{i = 1}^{N} \log \left(1 + \exp\left(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij}\right)\right) \tag{1.4} \end{eqnarray}\]
식 (1.4)을 각 회귀계수 \(\beta_0, \beta_1, \cdots, \beta_p\)에 대해 편미분하여 최적해를 얻는다. 이를 위해 주로 뉴턴-랩슨 알고리즘(Newton-Raphson algorithm)이나 quasi-Newton 알고리즘이 사용된다 (전치혁 2012).
회귀계수를 추정하는 함수를 fit_binary_logistic_regression()으로 구현하였다.
fit <- fit_binary_logistic_regression(student, y, c(x1, x2, x3), .reflevel = "보통")
print(fit)
#> $betas
#> (Intercept) x1 x2 x3
#> -30.521912 2.032088 3.471895 2.415513
#>
#> $hessian
#> (Intercept) x1 x2 x3
#> (Intercept) -1.6881640 -0.6963372 -12.650219 -2.595028
#> x1 -0.6963372 -0.6963372 -5.000818 -1.031367
#> x2 -12.6502186 -5.0008178 -96.025398 -17.990616
#> x3 -2.5950277 -1.0313675 -17.990616 -6.456938위 추정계수와 헤시안 행렬(Hessian matrix)를 이용하여 아래와 같이 결과를 요약할 수 있다.
fit_summary_df <- tibble(
term = names(fit$betas),
estimate = fit$betas,
std_error = sqrt((diag(solve(-fit$hessian)))),
statistic = estimate / std_error,
p_value = pnorm(if_else(statistic < 0, statistic, - statistic)) * 2,
odds_ratio = if_else(term == "(Intercept)", NA_real_, exp(estimate))
)| term | estimate | std_error | statistic | p_value | odds_ratio |
|---|---|---|---|---|---|
| (Intercept) | -30.52 | 18.03 | -1.69 | 0.09 | NA |
| x1 | 2.03 | 1.98 | 1.02 | 0.31 | 7.63 |
| x2 | 3.47 | 2.08 | 1.67 | 0.09 | 32.20 |
| x3 | 2.42 | 1.40 | 1.73 | 0.08 | 11.20 |
위 추정된 계수를 이용하여 사후확률을 계산하는 함수 회귀계수를 추정하는 함수를 posterior_binary_logistic_regression()로 구현하였다.
posterior_df <- posterior_binary_logistic_regression(
fit$betas, .new_data = student, c(x1, x2, x3),
.reflevel = "보통",
.poslevel = "우수"
)
prediction_df <- bind_cols(student, posterior_df)| id | x1 | x2 | x3 | y | .pred_우수 | .pred_보통 |
|---|---|---|---|---|---|---|
| 1 | 0 | 8 | 2 | 우수 | 0.89 | 0.11 |
| 2 | 1 | 7 | 1 | 우수 | 0.15 | 0.85 |
| 3 | 0 | 9 | 0 | 우수 | 0.67 | 0.33 |
| 4 | 1 | 6 | 4 | 우수 | 0.88 | 0.12 |
| 5 | 1 | 8 | 2 | 우수 | 0.98 | 0.02 |
| 6 | 0 | 7 | 3 | 우수 | 0.74 | 0.26 |
| 7 | 0 | 7 | 0 | 보통 | 0.00 | 1.00 |
| 8 | 1 | 6 | 1 | 보통 | 0.01 | 0.99 |
| 9 | 0 | 7 | 2 | 보통 | 0.20 | 0.80 |
| 10 | 0 | 8 | 1 | 보통 | 0.42 | 0.58 |
| 11 | 0 | 5 | 2 | 보통 | 0.00 | 1.00 |
| 12 | 1 | 8 | 0 | 보통 | 0.33 | 0.67 |
| 13 | 0 | 6 | 3 | 보통 | 0.08 | 0.92 |
| 14 | 1 | 7 | 2 | 보통 | 0.66 | 0.34 |
| 15 | 0 | 6 | 1 | 보통 | 0.00 | 1.00 |
종속변수가 셋 이상의 범주를 갖고 있으나 자연스러운 순서가 없는 경우, 기준범주 로짓모형이 널리 사용된다.
\(N\)개의 객체로 이루어진 학습데이터 \(\{(\mathbf{x}_i, y_i)\}_{i = 1, \cdots, N}\)를 아래와 같이 정의하자.
각 객체 \(i\)가 각 범주에 해당할 확률을 \(\pi_{ij}\)라 하자.
\[\begin{equation*} \pi_{ij} = P(y_i = j \, | \, \mathbf{x}_i), \, j = 1, \cdots, J \end{equation*}\]
이 때, 모든 \(i\)에 대하여
\[\begin{equation*} \sum_{j = 1}^{J} \pi_{ij} = 1 \end{equation*}\]
이 성립한다. 여기에서 범주 1을 기준 범주(reference category 혹은 baseline category)로 간주하여 범주별로 다음과 같은 회귀모형을 정의한다 (교재 (전치혁 2012)에는 범주 \(J\)를 기준 범주로 간주).
\[\begin{equation*} \log \left( \frac{\pi_{ij}}{\pi_{i1}} \right) = \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i, \, j = 2, \cdots, J \end{equation*}\]
이를 \(\pi_{ij}\)에 대해 풀면, 아래와 같은 해가 얻어진다 (Czepiel 2002).
\[\begin{equation} \begin{split} \pi_{ij} &= \frac{\exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)}{1 + \sum_{j = 2}^{J} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)}, \, j = 2, \cdots, J\\ \pi_{i1} &= \frac{1}{1 + \sum_{j = 2}^{J} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)} \end{split} \tag{2.1} \end{equation}\]
data(defecttype, package = "dmtr")| id | x1 | x2 | y |
|---|---|---|---|
| 1 | 0.09 | 5.02 | 1 |
| 2 | 0.10 | 5.01 | 1 |
| 3 | 0.12 | 4.94 | 1 |
| 4 | 0.12 | 5.12 | 1 |
| 5 | 0.12 | 5.03 | 1 |
| 6 | 0.12 | 4.94 | 2 |
| 7 | 0.10 | 5.13 | 2 |
| 8 | 0.10 | 4.87 | 1 |
| 9 | 0.10 | 5.13 | 2 |
| 10 | 0.11 | 4.94 | 3 |
| 11 | 0.11 | 4.93 | 3 |
| 12 | 0.09 | 5.02 | 3 |
| 13 | 0.10 | 5.01 | 3 |
| 14 | 0.09 | 4.94 | 3 |
| 15 | 0.10 | 5.12 | 2 |
| 16 | 0.12 | 4.93 | 2 |
| 17 | 0.10 | 5.00 | 1 |
| 18 | 0.09 | 5.01 | 3 |
위 모수 추정을 위해 최우추정법을 사용해보자. 우선, 종속변수를 변환한 지시변수를 아래와 같이 정의한다.
\[\begin{equation*} v_{ij} = \begin{cases} 1 & \text{ if } y_i = j\\ 0 & \text{ otherwise } \end{cases} \end{equation*}\]
이를 이용해 우도 함수를
\[\begin{eqnarray*} L &=& \prod_{i = 1}^{n} \prod_{j = 1}^{J} \left( \pi_{ij} \right)^{v_{ij}} \\ &=& \prod_{i = 1}^{n} \pi_{i1}^{1 - \sum_{j = 2}^{J} v_{ij}} \prod_{j = 2}^{J} \left( \pi_{ij} \right)^{v_{ij}}\\ &=& \prod_{i = 1}^{n} \frac{\pi_{i1}}{\pi_{i1}^{\sum_{j = 2}^{J} v_{ij}}} \prod_{j = 2}^{J} \left( \pi_{ij} \right)^{v_{ij}}\\ &=& \prod_{i = 1}^{n} \frac{\pi_{i1}}{\prod_{j = 2}^{J} \pi_{i1}^{v_{ij}}} \prod_{j = 2}^{J} \left( \pi_{ij} \right)^{v_{ij}}\\ &=& \prod_{i = 1}^{n} \pi_{i1} \prod_{j = 2}^{J} \left( \frac{\pi_{ij}}{\pi_{i1}} \right)^{v_{ij}} \end{eqnarray*}\]
와 같이 표현할 수 있으며, 여기에 식 (2.1)을 이용하면 아래와 같이 정리할 수 있다.
\[\begin{eqnarray*} L &=& \prod_{i = 1}^{n} \frac{1}{1 + \sum_{j = 2}^{J} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)} \prod_{j = 2}^{J} \left( \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right) \right)^{v_{ij}}\\ &=& \prod_{i = 1}^{n} \left( 1 + \sum_{j = 2}^{J} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right) \right)^{-1} \prod_{j = 2}^{J} \left( \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right) \right)^{v_{ij}} \end{eqnarray*}\]
이에 따라 로그 우도함수는 다음과 같이 정의된다.
\[\begin{equation} \log L = \sum_{i = 1}^{n} \left( - \log \left( 1 + \sum_{j = 2}^{J} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right) \right) + \sum_{j = 2}^{J} v_{ij} \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right) \right) \tag{2.2} \end{equation}\]
식 (2.2)을 각 계수에 대해 미분하면 아래와 같이 정리된다.
\[\begin{equation} \begin{split} \frac{\partial \log L}{\partial \beta_{0,j}} &= \sum_{i = 1}^{n} v_{ij} - \frac{\exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)}{1 + \sum_{j = 2}^{J} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)}\\ \frac{\partial \log L}{\partial \beta_{k,j}} &= \sum_{i = 1}^{n} v_{ij} x_{ik} - \frac{x_{ik} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)}{1 + \sum_{j = 2}^{J} \exp \left( \beta_{0,j} + \boldsymbol\beta_{j}^\top \mathbf{x}_i \right)}, \, k = 1, \cdots, p \end{split} \tag{2.3} \end{equation}\]
따라서, 명목 로지스틱 회귀분석은 식 (2.3)이 표현하는 \((J - 1) \times (p + 1)\)개의 미분식을 모두 0으로 만드는 계수값을 찾는 문제가 된다. 이에 대한 closed form solution은 존재하지 않으므로, 각종 알고리즘을 이용하여 해를 찾아야 한다. Newton-Raphson method에 의해 해를 찾는 방법은 Czepiel (2002) 에 보다 자세하게 설명되어 있다.
Table 2.1의 학습데이터에 대해 명목 로지스틱 회귀모형을 학습하여 범주를 추정한 결과는 아래와 같다.
회귀계수를 추정하는 함수를 fit_multinom_logistic_regression()으로 구현하였다.
fit <- fit_multinom_logistic_regression(defecttype, y, c(x1, x2), .reflevel = 3L,
.control = list(fnscale = -1, reltol = 1e-10, maxit = 100000))
print(fit)
#> $betas
#> 1 2
#> (Intercept) -63.65737 -117.04357
#> x1 101.59702 131.91194
#> x2 10.70778 20.61947
#>
#> $hessian
#> (Intercept):1 x1:1 x2:1 (Intercept):2 x1:2
#> (Intercept):1 -3.9801283 -0.41901950 -19.925885 1.9930964 0.21765391
#> x1:1 -0.4190195 -0.04459339 -2.096977 0.2176539 0.02397958
#> x2:1 -19.9258845 -2.09697702 -99.778149 10.0334423 1.09496537
#> (Intercept):2 1.9930964 0.21765392 10.033442 -2.8642989 -0.30647713
#> x1:2 0.2176539 0.02397958 1.094965 -0.3064771 -0.03312263
#> x2:2 10.0334492 1.09496616 50.522223 -14.3960085 -1.53952911
#> x2:2
#> (Intercept):1 10.033449
#> x1:1 1.094966
#> x2:1 50.522223
#> (Intercept):2 -14.396008
#> x1:2 -1.539529
#> x2:2 -72.371668종속변수가 3개 이상의 범주를 가지며, 각 범주 간에 서열이 있는 경우에 대한 로지스틱 회귀모형을 소개한다.
data(telconnection, package = "dmtr")| N | L | y |
|---|---|---|
| 25 | 5 | 3 |
| 25 | 10 | 3 |
| 25 | 20 | 2 |
| 25 | 30 | 1 |
| 32 | 5 | 3 |
| 32 | 10 | 3 |
| 32 | 20 | 2 |
| 32 | 30 | 1 |
| 42 | 5 | 1 |
| 42 | 10 | 3 |
| 42 | 20 | 1 |
| 42 | 30 | 1 |
객체 \(i\)가 범주 \(j\) 이하에 속할 누적확률을 \(\kappa_{ij}\)라 하자.
\[\begin{equation*} \kappa_{ij} = P(y_i \leq j \, | \, \mathbf{x}_i), \, j = 1, \cdots, J \end{equation*}\]
누적 로짓모형은 범주 누적확률의 로짓변환에 대한 선형 회귀모형이다.
\[\begin{equation} \log \left( \frac{\kappa_{ij}}{1 - \kappa_{ij}} \right) = \beta_{0,j} + \boldsymbol\beta^\top \mathbf{x}_i, \, j = 1, \cdots, J - 1 \tag{3.1} \end{equation}\]
식 (3.1)은 독립변수에 대한 계수 \(\boldsymbol\beta\)가 모든 범주에 대해 동일하며 절편(intercept) \(\beta_{0,j}\)만 범주에 따라 다른 비례 승산 모형(proportional odds model)이다. 즉, 범주에 관계없이 각 독립변수가 한 단위 증가할 때마다 로그 승산비는 동일하게 증가한다.
모형의 추정은 2절과 유사하게 다항분포를 사용한 최우추정법을 사용할 수 있다. 각 객체 \(i\)가 범주 \(j\)에 속할 확률은 아래와 같다.
\[\begin{equation*} \begin{split} \pi_{ij} &= \kappa_{ij} - \kappa_{i,j-1}\\ &= \frac{\exp (\beta_{0,j} + \boldsymbol\beta^\top \mathbf{x}_i)}{1 + \exp (\beta_{0,j} + \boldsymbol\beta^\top \mathbf{x}_i)} - \frac{\exp (\beta_{0,j-1} + \boldsymbol\beta^\top \mathbf{x}_i)}{1 + \exp (\beta_{0,j-1} + \boldsymbol\beta^\top \mathbf{x}_i)}, \, j = 2, \cdots, J - 1\\ & \\ \pi_{i1} &= \kappa_{i1}\\ &= \frac{\exp (\beta_{0,1} + \boldsymbol\beta^\top \mathbf{x}_i)}{1 + \exp (\beta_{0,1} + \boldsymbol\beta^\top \mathbf{x}_i)}\\ & \\ \pi_{iJ} &= 1 - \kappa_{i,J-1}\\ &= 1 - \frac{\exp (\beta_{0,J-1} + \boldsymbol\beta^\top \mathbf{x}_i)}{1 + \exp (\beta_{0,J-1} + \boldsymbol\beta^\top \mathbf{x}_i)} \end{split} \tag{3.2} \end{equation*}\]
로그 우도함수는
\[\begin{equation*} \sum_{i = 1}^{n} \sum_{j = 1}^{J} y_i \log \pi_{ij} \end{equation*}\]
이며, 이에 위에서 정리한 \(\pi_{ij}\)식을 대입하여 전개할 수 있다. 이 로그 우도함수는 concave 함수이므로(Pratt 1981), 각 계수에 대해 편미분하여 0이 되도록 하는 값을 구하는 방식으로 회귀모형을 추정할 수 있다.
회귀계수를 추정하는 함수를 fit_ordinal_logistic_regression()으로 구현하였다.
fit <- fit_ordinal_logistic_regression(
telconnection, y, c(N, L)
)
print(fit)
#> $par
#> [1] -13.0270423 -11.3909924 0.2234882 0.2996985
#>
#> $value
#> [1] -6.441254
#>
#> $counts
#> function gradient
#> 449 NA
#>
#> $convergence
#> [1] 0
#>
#> $message
#> NULL
#>
#> $hessian
#> [,1] [,2] [,3] [,4]
#> [1,] -1.3670843 0.6006743 -24.21539 -15.81477
#> [2,] 0.6006743 -1.4690612 -27.55876 -12.06814
#> [3,] -24.2153900 -27.5587587 -1707.10486 -843.13495
#> [4,] -15.8147679 -12.0681364 -843.13495 -559.77999