In youngroklee-ml/dmtr: Functions and Data for "Data Mining Techniques"

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(dmtr)
library(dplyr)

이분 로지스틱 회귀모형 {#binary-logistic-reg-model}

이분 로지스틱 회귀모형은 종속변수가 2가지 범주를 취하는 경우에 사용된다.

$N$개의 객체로 이루어진 학습데이터 ${(\mathbf{x}i, y_i)}{i = 1, \cdots, N}$를 아래와 같이 정의하자.

• $\mathbf{x}i \in \mathbb{R}^p$: $p$개의 독립변수로 이루어진 벡터 ($\mathbf{x}_i = [x{i1} \, x_{i2} \, \cdots \, x_{ip}]^\top$)
• $y_i$: 0 혹은 1의 값을 갖는 이분형 지시변수 (indicator variable)

$\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)} (#eq:logistic-function) \end{eqnarray}

여기에서 $\boldsymbol\beta \in \mathbb{R}^{p}$는 $\mathbf{x}_i$와 동일한 차원의 벡터이다 ($\boldsymbol\beta = [\beta_1 \, \beta_2 \, \cdots \, \beta_p]^\top$).

식 \@ref(eq:logistic-function)는 모든 $\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 (#eq:logit-transform) \end{eqnarray}

식 \@ref(eq:logit-transform)에서 확률 $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} (#eq:binary-logistic-latent-interpret) \end{equation}

식 \@ref(eq:binary-logistic-latent-interpret)에서 $\varepsilon_i$는 표준로지스틱분포(standard logistic distribution)을 따른다.

데이터

data(student, package = "dmtr")

knitr::kable(
  student,
  align = c('r', 'r', 'r', 'r', 'r'),
  caption = '우수/보통 학생에 대한 설문조사 결과'
)

회귀계수 추정 {#binary-logistic-reg-estimation}

로지스틱 모형에서 회귀계수의 추정을 위해서 주로 최우추정법(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) (#eq:binary-logistic-reg-loglik) \end{eqnarray}

식 \@ref(eq:binary-logistic-reg-loglik)을 각 회귀계수 $\beta_0, \beta_1, \cdots, \beta_p$에 대해 편미분하여 최적해를 얻는다. 이를 위해 주로 뉴턴-랩슨 알고리즘(Newton-Raphson algorithm)이나 quasi-Newton 알고리즘이 사용된다 [@jun2012datamining].

회귀계수를 추정하는 함수를 fit_binary_logistic_regression()으로 구현하였다.

fit <- fit_binary_logistic_regression(student, y, c(x1, x2, x3), .reflevel = "보통")
print(fit)

위 추정계수와 헤시안 행렬(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))
)

knitr::kable(
  fit_summary_df,
  booktabs = TRUE,
  digits = 2L,
  caption = "우수/보통 학생 설문조사 데이터에 대한 Logistic Regression 결과"
)

위 추정된 계수를 이용하여 사후확률을 계산하는 함수 회귀계수를 추정하는 함수를 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)

knitr::kable(
  prediction_df,
  booktabs = TRUE,
  align = rep('r', ncol(prediction_df)),
  digits = 2L,
  caption = "우수/보통 학생 설문조사 데이터에 대한 범주 추정"
)

명목 로지스틱 회귀모형 {#nominal-logistic-regression}

종속변수가 셋 이상의 범주를 갖고 있으나 자연스러운 순서가 없는 경우, 기준범주 로짓모형이 널리 사용된다.

$N$개의 객체로 이루어진 학습데이터 ${(\mathbf{x}i, y_i)}{i = 1, \cdots, N}$를 아래와 같이 정의하자.

• $\mathbf{x}i \in \mathbb{R}^p$: $p$개의 독립변수로 이루어진 벡터 ($\mathbf{x}_i = [x{i1} \, x_{i2} \, \cdots \, x_{ip}]^\top$)
• $J$: 범주 수
• $y_i$: 객체 $i$에 대한 종속변수값 $\in {1, 2, \cdots, J}$

각 객체 $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)로 간주하여 범주별로 다음과 같은 회귀모형을 정의한다 (교재 [@jun2012datamining]에는 범주 $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}$에 대해 풀면, 아래와 같은 해가 얻어진다 [@czepiel2002maximum].

\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} (#eq:multi-nominal-prob-sol) \end{equation}

데이터

data(defecttype, package = "dmtr")

knitr::kable(defecttype, booktabs = TRUE,
  align = c('r', 'r', 'r', 'r'),
  caption = '공정변수-불량 종류 데이터')

회귀계수 추정 {#multinom-logistic-reg-estimation}

위 모수 추정을 위해 최우추정법을 사용해보자. 우선, 종속변수를 변환한 지시변수를 아래와 같이 정의한다.

\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}

와 같이 표현할 수 있으며, 여기에 식 \@ref(eq:multi-nominal-prob-sol)을 이용하면 아래와 같이 정리할 수 있다.

\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) (#eq:multi-nominal-logit-loglik) \end{equation}

식 \@ref(eq:multi-nominal-logit-loglik)을 각 계수에 대해 미분하면 아래와 같이 정리된다.

\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} (#eq:multi-nominal-logit-loglik-diff) \end{equation}

따라서, 명목 로지스틱 회귀분석은 식 \@ref(eq:multi-nominal-logit-loglik-diff)이 표현하는 $(J - 1) \times (p + 1)$개의 미분식을 모두 0으로 만드는 계수값을 찾는 문제가 된다. 이에 대한 closed form solution은 존재하지 않으므로, 각종 알고리즘을 이용하여 해를 찾아야 한다. Newton-Raphson method에 의해 해를 찾는 방법은 @czepiel2002maximum 에 보다 자세하게 설명되어 있다.

Table \@ref(tab:nominal-logistic-reg-train-data-print)의 학습데이터에 대해 명목 로지스틱 회귀모형을 학습하여 범주를 추정한 결과는 아래와 같다.

회귀계수를 추정하는 함수를 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)

서열 로지스틱 회귀모형 {#ordinal-logistic-regression}

종속변수가 3개 이상의 범주를 가지며, 각 범주 간에 서열이 있는 경우에 대한 로지스틱 회귀모형을 소개한다.

데이터

data(telconnection, package = "dmtr")

knitr::kable(telconnection, booktabs = TRUE,
  align = c('r', 'r', 'r'),
  label = c('잡음(N)', '손실(L)', '만족도(y)'),
  caption = '성능변수에 따른 통신 만족도')

누적 로짓모형 {#cumulative-logit-model}

객체 $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 (#eq:cumulative-logit) \end{equation}

식 \@ref(eq:cumulative-logit)은 독립변수에 대한 계수 $\boldsymbol\beta$가 모든 범주에 대해 동일하며 절편(intercept) $\beta_{0,j}$만 범주에 따라 다른 비례 승산 모형(proportional odds model)이다. 즉, 범주에 관계없이 각 독립변수가 한 단위 증가할 때마다 로그 승산비는 동일하게 증가한다.

모형의 추정은 \@ref(nominal-logistic-regression)절과 유사하게 다항분포를 사용한 최우추정법을 사용할 수 있다. 각 객체 $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} (#eq:cumulative-logit-prob) \end{equation}

로그 우도함수는

\begin{equation} \sum_{i = 1}^{n} \sum_{j = 1}^{J} y_i \log \pi_{ij} \end{equation}

이며, 이에 위에서 정리한 $\pi_{ij}$식을 대입하여 전개할 수 있다. 이 로그 우도함수는 concave 함수이므로[@pratt1981concavity], 각 계수에 대해 편미분하여 0이 되도록 하는 값을 구하는 방식으로 회귀모형을 추정할 수 있다.

회귀계수를 추정하는 함수를 fit_ordinal_logistic_regression()으로 구현하였다.

fit <- fit_ordinal_logistic_regression(
  telconnection, y, c(N, L)
)
print(fit)

youngroklee-ml/dmtr documentation built on June 12, 2022, 6:24 p.m.