Chapter 6 로지스틱 회귀분석

로지스틱 회귀분석(logistic regression)은 종속변수가 통상 2개의 범주(있음/없음, 불량/양호, 합격/불합격 등)를 다루는 모형을 지칭하나, 3개 이상의 범주를 다루기도 한다. 후자의 경우는 다시 서열형(ordinal) 데이터와 명목형(nominal) 데이터인 경우에 따라 서로 다른 모형이 사용된다.

6.1 필요 R 패키지 설치

본 장에서 필요한 R 패키지들은 아래와 같다.

package	version
tidyverse	1.3.1
stats	4.1.3
nnet	7.3-16
MASS	7.3-54
VGAM	1.1-5

6.2 이분 로지스틱 회귀모형

6.2.1 기본 R 스크립트

train_df <- tribble(
  ~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, "보통"
) %>%
  mutate(y = factor(y, levels = c("보통", "우수")))

knitr::kable(train_df, booktabs = TRUE,
             align = c('r', 'r', 'r', 'r', 'r'),
             col.names = c('객체번호', '아침식사여부($x_1$)', '수면시간($x_2$)', '서클활동시간($x_3$)', '범주(y)'),
             caption = '우수/보통 학생에 대한 설문조사 결과')

Table 6.1: 우수/보통 학생에 대한 설문조사 결과

객체번호	아침식사여부(\(x_1\))	수면시간(\(x_2\))	서클활동시간(\(x_3\))	범주(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	보통

Table 6.1와 같이 세 개의 독립변수 \(x_1\), \(x_2\), \(x_3\)와 이분형 종속변수 \(y\)의 관측값(보통 = 0, 우수 = 1)으로 이루어진 15개의 학습표본을 train_df라는 data frame에 저장한다.

아래와 같이 glm 함수를 이용하여 로지스틱 회귀모형을 간편하게 추정할 수 있다.

glm_fit <- glm(y ~ x1 + x2 + x3, family = binomial(link = "logit"), data = train_df)

knitr::kable(
  glm_fit %>%  broom::tidy(),
  booktabs = TRUE,
  caption = "위 우수/보통 학생 설문조사 데이터에 대한 Logistic Regression 결과"
#  caption = "Table \\@ref(tab:binary-logistic-reg-train-data)에 대한 Logistic Regression 결과"
  )

Table 6.2: 위 우수/보통 학생 설문조사 데이터에 대한 Logistic Regression 결과

term	estimate	std.error	statistic	p.value
(Intercept)	-30.510836	18.018256	-1.693329	0.0903929
x1	2.031278	1.983692	1.023989	0.3058406
x2	3.470671	2.074978	1.672631	0.0944000
x3	2.414387	1.396372	1.729043	0.0838015

Table 6.2은 추정된 회귀계수 추정치 estmate과 그 표준오차 std.error, 표준화(standardized)된 회귀계수값 statistic (= estmate / std.error), 그리고 귀무가설 \(H_0\): statistic = 0 에 대한 유의확률 p.value를 보여준다.

6.2.2 회귀모형

이분 로지스틱 회귀모형은 종속변수가 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)} \tag{6.1} \end{eqnarray}\]

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

식 (6.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{6.2} \end{eqnarray}\]

식 (6.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{6.3} \end{equation}\]

식 (6.3)에서 \(\varepsilon_i\)는 표준로지스틱분포(standard logistic distribution)을 따른다.

6.2.3 회귀계수 추정

로지스틱 모형에서 회귀계수의 추정을 위해서 주로 최우추정법(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{6.4} \end{eqnarray}\]

식 (6.4)을 각 회귀계수 \(\beta_0, \beta_1, \cdots, \beta_p\)에 대해 편미분하여 최적해를 얻는다. 이를 위해 주로 뉴턴-랩슨 알고리즘(Newton-Raphson algorithm)이나 quasi-Newton 알고리즘이 사용되나 (전치혁 2012), 본 장에서는 우선 안정성은 떨어지지만 보다 간편한 방법으로 경사하강법(gradient descent)을 소개한다.

6.2.3.1 경사하강법

식 (6.1)과 \(P(y_i = 0 \,|\, \mathbf{x}_i) = 1 - P_i\), 그리고

\[\begin{equation*} \frac{\exp(z)}{1 + \exp(z)} = \frac{1}{1 + \exp(-z)} \end{equation*}\]

임을 고려하면 아래와 같이 범주확률모형을 정의할 수 있다.

\[\begin{equation*} P(y = y_i \,|\, \mathbf{x}_i, \beta_0, \boldsymbol\beta) = \frac{1}{1 + \exp\left((1 - 2y_i)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)} \end{equation*}\]

이에 따라 로그우도함수 (6.4)는 아래와 같이 정리된다.

\[\begin{equation*} \log \prod_{i = 1}^{N} P(y = y_i \,|\, \mathbf{x}_i, \beta_0, \boldsymbol\beta) = - \sum_{i = 1}^{N} \log \left(1 + \exp\left((1 - 2y_i)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)\right) \end{equation*}\]

위 로그우도함수를 최대화하는 문제는 아래 함수를 최소화하는 문제와 동일하다.

\[\begin{equation} f(\beta_0, \boldsymbol\beta) = \sum_{i = 1}^{N} \log \left(1 + \exp\left((1 - 2y_i)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)\right) \tag{6.5} \end{equation}\]

경사하강법에 따라 아래의 과정을 통해 회귀계수를 추정할 수 있다.

1. 임의의 값으로 \(\beta_0, \beta_1, \cdots, \beta_j\)의 초기 추정값을 설정한다.
2. 식 (6.5)을 각 회귀변수에 대해 편미분한 미분값을 구한다.
3. 2의 값에 학습률(step size)을 곱한 만큼 회귀계수 추정값을 이동시킨다. 방향은 미분값의 반대방향.
4. 수렴할 때까지 2-3의 과정을 반복한다.

여기에서 식 (6.5)의 각 회귀변수에 대한 편미분식은 아래와 같다.

\[\begin{eqnarray*} \frac{\partial f}{\partial \beta_0} &=& \sum_{i = 1}^{N} (1 - 2y_i) \frac{\exp\left((1 - 2y_i)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)}{1 + \exp\left((1 - 2y_i)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)}\\ &=& \sum_{i = 1}^{N} \frac{1 - 2y_i}{1 + \exp\left((2y_i - 1)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)} \end{eqnarray*}\]

\[\begin{eqnarray*} \frac{\partial f}{\partial \beta_j} &=& \sum_{i = 1}^{N} (1 - 2y_i)x_{ij} \frac{\exp\left((1 - 2y_i)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)}{1 + \exp\left((1 - 2y_i)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)}\\ &=& \sum_{i = 1}^{N} \frac{(1 - 2y_i)x_{ij}}{1 + \exp\left((2y_i - 1)(\beta_0 + \sum_{j = 1}^{p} \beta_j x_{ij})\right)} \end{eqnarray*}\]

따라서, 회귀계수 추정값을 이동시키는 함수 update_beta를 아래와 같이 구현할 수 있다.

update_beta <- function(x, y, beta0 = 0, beta = rep(0, dim(x)[2]), alpha = 0.01) {
  # 변미분식의 분모
  denominator <- 1 + exp((2 * y - 1) * (beta0 + (x %*% beta)))

  # intercept 이동량 계산
  beta0_numerator <- 1 - 2 * y
  beta0_update = sum(beta0_numerator / denominator)
  
  # intercept 외 회귀계수 이동량 계산
  beta_numerator <- sweep(x, MARGIN = 1, STATS = 1 - 2 * y, FUN = "*")
  beta_update = apply(beta_numerator, MARGIN = 2, 
                      function(x) sum(x / denominator))

  # 회귀계수 이동
  beta0 <- beta0 - alpha * beta0_update
  beta <- beta - alpha * beta_update

  return(list(beta0 = beta0, beta = beta))
}

위의 함수를 이용하여 아래 estimate_beta처엄 수렴할 때까지 회귀계수 추정값을 계속 이동시킨다. 본 경사하강법은 학습률 파라미터 alpha값에 따라 민감한 단점이 있으며, 특히 alpha값을 크게 설정할 경우에는 추정값이 수렴하지 않고 오히려 실제값에서 계속 멀어지는 현상이 발생하기도 한다. 이러한 단점을 보완하기 위한 여러 방법이 있으나, 본 장에서 자세한 설명은 생략하기로 한다.

caltculate_loglik <- function(x, y, 
                              beta0 = 0, 
                              beta = rep(0, dim(x)[2])) {
  sum(y * (beta0 + (x %*% beta))) - 
    sum(log(1 + exp(beta0 + (x %*% beta))))
}

estimate_beta <- function(x, y, 
                          beta0 = 0, 
                          beta = rep(0, dim(x)[2]), 
                          alpha = 0.01, 
                          conv_threshold = 1e-5, 
                          max_iter = 1e+5) {
  new_beta0 <- beta0
  new_beta <- beta
  conv <- FALSE
  
  i_iter <- 0
  while(i_iter < max_iter) {
    res <- update_beta(x, y, new_beta0, new_beta, alpha)
    
    if(abs(caltculate_loglik(x, y, beta0, beta)
           - caltculate_loglik(x, y, res$beta0, res$beta))
       < conv_threshold) {
      conv <- TRUE
      break
      }
    
    new_beta0 <- res$beta0
    new_beta <- res$beta
    
    i_iter <- i_iter + 1
  }
  
  return(list(conv = conv, beta0 = new_beta0, beta = new_beta))
}

위에서 정의한 함수를 이용하여 Table 6.1의 학습표본에 대한 로지스틱 회귀모형을 추정해보자.

res <- estimate_beta(train_df[, c("x1", "x2", "x3")] %>% as.matrix(), 
                     train_df$y %>% as.numeric() - 1,
                     alpha = 0.015,
                     conv_threshold = 1e-6)

print(res)

## $conv
## [1] FALSE
## 
## $beta0
## [1] -30.39189
## 
## $beta
##       x1       x2       x3 
## 2.022808 3.457019 2.405969

위 회귀계수 추정값은 R 함수 glm을 이용한 추정값(Table 6.2)과 유사함을 볼 수 있다.

6.2.3.2 반복재가중최소제곱법

R의 glm 함수는 반복재가중최소제곱법(iteratively rewighted least squares; IRLS 혹은 IWLS)을 사용한다. 이는 선형회귀식

\[\begin{equation*} logit(y_i) = \beta_0 + \boldsymbol\beta^\top \mathbf{x}_i + \varepsilon_i \end{equation*}\]

을 추정하는 방법인데, 여기에서 \(y_i\)는 0 혹은 1이므로, \(logit(y_i)\)는 \(-\infty\) 혹은 \(\infty\)가 되어 회귀식을 추정할 수 없다. 따라서, 식 (6.1)에 설명된 로지스틱 함수

\[\begin{equation*} P = \frac{\exp(\beta_0 + \boldsymbol\beta^\top \mathbf{x})}{1 + \exp(\beta_0 + \boldsymbol\beta^\top \mathbf{x})} \end{equation*}\]

와 테일러 급수(Taylor series)를 이용하여 \(logit(y)\)에 대한 근사함수를 아래와 같이 얻는다.

\[\begin{eqnarray*} g(y) &=& logit(P) + (y - P) \frac{\partial logit(P)}{\partial P}\\ &=& \log \frac{P}{1 - P} + (y - P) \left( \frac{1}{P} + \frac{1}{1 - P} \right) \end{eqnarray*}\]

그리고 아래 선형회귀식을 추정한다.

\[\begin{equation*} g(y_i) = \beta_0 + \boldsymbol\beta^\top \mathbf{x}_i + \varepsilon_i \end{equation*}\]

여기에서 오차항 \(\varepsilon_i\)의 분산은 추정된 확률 \(P_i\)에 따라 다르므로, 통상적 최소자승법(ordinary least squares; OLS) 대신 오차항의 분산이 동일해지도록 객체마다 가중치를 부여하는 가중최소자승법(weighted least squares; WLS)을 사용한다. 로지스틱 회귀모형에서 각 객체의 가중치는

\[\begin{equation*} w_i = P_i (1 - P_i) \end{equation*}\]

가중치와 회귀계수 추정값은 상호 영향을 미치므로, 수렴할 때까지 반복적으로 가중치와 회귀계수 추정값을 변화시키면서 최종 추정값을 찾아가는 방법이다.

우선 회귀계수 추정값이 주어졌을 때 각 객체에 대한 확률값 \(P_i\)와 가중치 \(w_i\)를 구하는 함수 calculate_weight를 아래와 같이 구현해보자.

calculate_weight <- function(x, beta0 = 0, 
                             beta = rep(0, dim(x)[2])) {
  # 각 객체의 y값이 1일 확률
  P <- (1 + exp(- beta0 - (x %*% beta)))^(-1) %>% drop()

  # 가중치 계산
  w <- P * (1 - P)

  return(list(P = P, w = w))
}

그리고 확률추정값과 가중치가 주어졌을 때 회귀계수를 구하는 함수 calculate_beta를 아래와 같이 구현해보자. 여기서 회귀계수를 구하는 부분은 R의 선형회귀분석함수 lm을 사용한다.

calculate_beta <- function(x, y, P, w) {
  # 추정확률값이 0 이나 1인 경우 
  # 여전히 logit 함수가 정의되지 않으므로 회귀계수 결정에서 제외
  logit_derivative <- 1/P + 1/(1 - P)
  is_good <- !is.nan(logit_derivative)
  
  # 모든 객체에 대한 추정확률이 0 이나 1인 경우 회귀계수 추정 불가능
  if(all(!is_good)) return(NULL)
  
  # 테일러 급수 계산
  g_y <- log(P[is_good]) - 
    log(1 - P[is_good]) + 
    (y[is_good] - P[is_good]) * logit_derivative
  
  # 가중치최소자승법을 이용한 추정
  df <- bind_cols(as_tibble(x) %>% 
                    `colnames<-`(colnames(x)), 
                  tibble(g_y = g_y))
  lm(g_y ~ ., data = df, subset = is_good, weights = w)
}

위에서 정의한 두 함수 calculate_weight과 calculate_beta를 반복적으로 사용하여 Table 6.1의 학습표본에 대한 로지스틱 회귀모형을 추정해보자. 모든 객체의 가중치 변화량이 1/10000 보다 작을 경우 모형추정이 수렴한 것으로 간주하도록 하자.

X <- train_df[, c("x1", "x2", "x3")] %>% as.matrix()
y <- train_df$y %>% as.numeric() - 1

weight <- calculate_weight(X)
for(i in 1:10) {
  wls_fit <- calculate_beta(X, y, weight$P, weight$w)
  
  if(is.null(wls_fit)) {break}
  
  new_weight <- calculate_weight(x = X, 
                                 beta0 = coef(wls_fit)[1], 
                                 beta = coef(wls_fit)[-1])
  
  if(max(abs(new_weight$w - weight$w)) < 1e-4) {break}
  
  weight <- new_weight
}

coef(wls_fit)

## (Intercept)          x1          x2          x3 
##  -30.510837    2.031278    3.470671    2.414387

위 스크립트를 실행시킨 결과 7번째 반복수행에서 결과가 수렴하였으며, 해당 결과는 glm 함수를 사용하였을 때의 결과 (Table 6.2)과 매우 근사함을 확인할 수 있다.

6.3 명목 로지스틱 회귀모형

6.3.1 기본 R 스크립트

train_df <- tribble(
  ~id, ~x1, ~x2, ~y,
  1, 0.09, 5.02, 1,
  2, 0.1, 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.1, 5.13, 2,
  8, 0.1, 4.87, 1,
  9, 0.1, 5.13, 2,
  10, 0.11, 4.94, 3,
  11, 0.11, 4.93, 3,
  12, 0.09, 5.02, 3,
  13, 0.1, 5.01, 3,
  14, 0.09, 4.94, 3,
  15, 0.1, 5.12, 2,
  16, 0.12, 4.93, 2,
  17, 0.1, 5, 1,
  18, 0.09, 5.01, 3
) %>%
  mutate(y = factor(y, levels = c(1, 2, 3)))

knitr::kable(train_df, booktabs = TRUE,
             align = c('r', 'r', 'r', 'r'),
             col.names = c('객체번호', '$x_1$', '$x_2$', '불량범주($y$)'),
             caption = '공정변수-불량 종류 데이터')

Table 6.3: 공정변수-불량 종류 데이터

객체번호	\(x_1\)	\(x_2\)	불량범주(\(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

Table 6.3와 같이 두 개의 독립변수 \(x_1\), \(x_2\)에 따라 세 종류의 불량 (\(y = 1, 2, 3\))$이 발생함을 알았다면, 아래와 같이 nnet 패키지의 multinom 함수를 이용하여 공정변수에 따른 불량 종류를 분류하기 위한 로지스틱 회귀모형을 간편하게 추정할 수 있다.

multinom_fit <- nnet::multinom(
  y ~ x1 + x2, 
  data = train_df, 
  maxit = 1000)

## # weights:  12 (6 variable)
## initial  value 19.775021 
## iter  10 value 17.825831
## iter  20 value 16.489924
## iter  30 value 16.116751
## iter  40 value 16.044223
## iter  50 value 16.025259
## iter  60 value 16.024629
## iter  70 value 16.016395
## final  value 16.015185 
## converged

knitr::kable(
  multinom_fit %>%  broom::tidy(exponentiate = FALSE),
  booktabs = TRUE,
  caption = '위 공정변수-불량종류 데이터에 대한 Logistic Regression 결과'
#  caption = 'Table \\@ref(tab:nominal-logistic-reg-train-data)에 대한 Logistic Regression 결과'
  )

Table 6.4: 위 공정변수-불량종류 데이터에 대한 Logistic Regression 결과

y.level	term	estimate	std.error	statistic	p.value
2	(Intercept)	-53.924661	45.402992	-1.1876896	0.2349557
2	x1	31.545749	62.162158	0.5074751	0.6118215
2	x2	9.992311	8.479305	1.1784352	0.2386231
3	(Intercept)	64.191145	57.282533	1.1206059	0.2624556
3	x1	-101.817613	65.516143	-1.5540844	0.1201643
3	x2	-10.810865	10.915881	-0.9903795	0.3219887

multinom 함수는 범주 y의 값 1, 2, 3중 첫번째 값인 1을 기준범주(reference category)로 사용한다.

6.3.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\))
• \(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)로 간주하여 범주별로 다음과 같은 회귀모형을 정의한다 (교재 (전치혁 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{6.6} \end{equation}\]

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

\[\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*}\]

와 같이 표현할 수 있으며, 여기에 식 (6.6)을 이용하면 아래와 같이 정리할 수 있다.

\[\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{6.7} \end{equation}\]

식 (6.7)을 각 계수에 대해 미분하면 아래와 같이 정리된다.

\[\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{6.8} \end{equation}\]

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

Table 6.3의 학습데이터에 대해 명목 로지스틱 회귀모형을 학습하여 범주를 추정한 결과는 아래와 같다.

multinom_fit <- nnet::multinom(
  y ~ x1 + x2, 
  data = train_df, 
  maxit = 1000)

## # weights:  12 (6 variable)
## initial  value 19.775021 
## iter  10 value 17.825831
## iter  20 value 16.489924
## iter  30 value 16.116751
## iter  40 value 16.044223
## iter  50 value 16.025259
## iter  60 value 16.024629
## iter  70 value 16.016395
## final  value 16.015185 
## converged

predict_df <- predict(multinom_fit, train_df, type = "probs") %>% 
  as_data_frame() %>%
  `colnames<-`(c("p1", "p2", "p3")) %>%
  mutate(pred_class = predict(multinom_fit, train_df, type = "class"))

## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.

bind_cols(train_df, predict_df) %>%
  select(-x1, -x2) %>%
  knitr::kable(
    booktabs = TRUE,
    align = c('r', 'r', 'r', 'r', 'r', 'r'),
    col.names = c('객체번호', '불량범주 $y_i$', 
                  '$\\pi_{i1}$', '$\\pi_{i2}$', '$\\pi_{i3}$', '추정범주 $\\hat{y}_i$'),
    caption = '명목 로지스틱 회귀모형 범주 추정 결과',
    digits = 3
  )

Table 6.5: 명목 로지스틱 회귀모형 범주 추정 결과

객체번호	불량범주 \(y_i\)	\(\pi_{i1}\)	\(\pi_{i2}\)	\(\pi_{i3}\)	추정범주 \(\hat{y}_i\)
1	1	0.283	0.112	0.604	3
2	1	0.425	0.209	0.365	1
3	1	0.589	0.271	0.141	1
4	1	0.262	0.729	0.009	2
5	1	0.450	0.509	0.041	2
6	2	0.589	0.271	0.141	1
7	2	0.349	0.570	0.082	2
8	1	0.199	0.024	0.777	3
9	2	0.349	0.570	0.082	2
10	3	0.501	0.168	0.331	1
11	3	0.490	0.149	0.361	1
12	3	0.283	0.112	0.604	3
13	3	0.425	0.209	0.365	1
14	3	0.160	0.029	0.811	3
15	2	0.365	0.540	0.095	2
16	2	0.595	0.247	0.158	1
17	1	0.416	0.186	0.398	1
18	3	0.268	0.096	0.636	3

nnet 패키지 외에도 glmnet, mlogit, VGAM 등의 R 패키지들을 사용해 명목형 로지스틱 회귀모형을 추정할 수 있다.

VGAM::vglm(y ~ x1 + x2,
           data = train_df,
           family = VGAM::multinomial)

## 
## Call:
## VGAM::vglm(formula = y ~ x1 + x2, family = VGAM::multinomial, 
##     data = train_df)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2          x1:1          x1:2 
##     -64.56378    -118.15274     102.65063     133.29235 
##          x2:1          x2:2 
##      10.86829      20.81307 
## 
## Degrees of Freedom: 36 Total; 30 Residual
## Residual deviance: 32.03007 
## Log-likelihood: -16.01503 
## 
## This is a multinomial logit model with 3 levels

6.4 서열 로지스틱 회귀모형

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

6.4.1 기본 R 스크립트

train_df <- tribble(
  ~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
) %>%
  mutate(y = as.ordered(y))

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

Table 6.6: 성능변수에 따른 통신 만족도

잡음(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

Table 6.6은 벨 연구소에서 한 통신장치에 대하여 실시한 조사 결과를 나타낸 것이다. 주요 성능변수인 회선잡음(circuit noise: N)과 소리크기 손실(loudness loss: L)이 이용자의 주관적인 만족도에 미치는 영향을 분석하기 위한 것이다. 만족도는 원결과(Cavanaugh, Hatch, and Sullivan 1976)를 가공하여 다음과 같이 3가지로 분류하였다.

\[\begin{equation*} y = \begin{cases} 1 & \mbox{good}\\ 2 & \mbox{fair}\\ 3 & \mbox{poor} \end{cases} \end{equation*}\]

본 장에서는 두 가지 모형을 다룬다. 우선 누적 로짓모형(cumulative logit model)은 아래와 같이 MASS 패키지의 polr 함수를 사용하여 추정할 수 있다.

MASS::polr(y ~ N + L, data = train_df) %>%
  broom::tidy()

## 
## Re-fitting to get Hessian

## # A tibble: 4 x 5
##   term  estimate std.error statistic coef.type  
##   <chr>    <dbl>     <dbl>     <dbl> <chr>      
## 1 N       -0.224     0.146     -1.53 coefficient
## 2 L       -0.300     0.137     -2.19 coefficient
## 3 1|2    -13.0       6.46      -2.02 scale      
## 4 2|3    -11.4       6.17      -1.85 scale

인근범주 로짓모형(adjacent-categories logit model)은 아래와 같이 VGAM 패키지의 vglm 함수를 이용하여 추정할 수 있다.

VGAM::vglm(y ~ N + L,
           data = train_df,
           family = VGAM::acat(reverse = TRUE))

## 
## Call:
## VGAM::vglm(formula = y ~ N + L, family = VGAM::acat(reverse = TRUE), 
##     data = train_df)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2           N:1           N:2 
##  -12.94227208   -6.10597713    0.31431599    0.04643039 
##           L:1           L:2 
##    0.17500408    0.29578067 
## 
## Degrees of Freedom: 24 Total; 18 Residual
## Residual deviance: 11.67769 
## Log-likelihood: -5.838847 
## 
## This is an adjacent categories model with 3 levels

6.4.2 누적 로짓모형

객체 \(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{6.9} \end{equation}\]

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

모형의 추정은 6.3.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{6.10} \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이 되도록 하는 값을 구하는 방식으로 회귀모형을 추정할 수 있다.

polr_fit <- MASS::polr(y ~ N + L, data = train_df)

print(polr_fit)

## Call:
## MASS::polr(formula = y ~ N + L, data = train_df)
## 
## Coefficients:
##          N          L 
## -0.2236292 -0.2998833 
## 
## Intercepts:
##       1|2       2|3 
## -13.03527 -11.39902 
## 
## Residual Deviance: 12.8825 
## AIC: 20.8825

위와 같이 polr 함수 실행 시 얻어지는 각 변수들에 대한 계수들의 부호는 교재(전치혁 2012)의 내용과 반대인데, 이는 polr 함수는 아래와 같은 모형을 추정하기 때문이다.

\[\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 \end{equation*}\]

위 모형에서 polr 함수 실행 결과 추정된 절편값은 \(\beta_{0,1} = -13.0352721\), \(\beta_{0,2} = -11.3990207\) 이며, 두 변수 \(N\), \(L\)에 대한 회귀계수는 각각 -0.2236292, -0.2998833로 추정된다.

추정된 회귀계수를 식 (6.10)에 대입하면 각 객체 \(i\)가 각 범주 \(j\)에 속할 확률을 Table 6.7와 같이 얻을 수 있다. 아래 R 스크립트에서 사용한 predict라는 함수가 해당 계산을 수행한다.

predict_df <- predict(polr_fit, train_df, type = "probs") %>% 
  as_data_frame() %>%
  `colnames<-`(c("p1", "p2", "p3")) %>%
  mutate(pred_class = predict(polr_fit, train_df, type = "class"))

bind_cols(train_df, predict_df) %>%
  knitr::kable(
    booktabs = TRUE,
    align = c('r', 'r', 'r', 'r', 'r', 'r', 'r'),
    col.names = c('잡음(N)', '손실(L)', '실제범주 $y_i$', 
                  '$\\pi_{i1}$', '$\\pi_{i2}$', '$\\pi_{i3}$', '추정범주 $\\hat{y}_i$'),
    caption = '위 통신 만족도 데이터에 대한 누적 로짓모형의 추정범주',
#    caption = 'Table \\@ref(tab:ordinal-logistic-reg-train-data)에 대한 누적 로짓모형의 추정범주',
    digits = 4
  )

Table 6.7: 위 통신 만족도 데이터에 대한 누적 로짓모형의 추정범주

잡음(N)	손실(L)	실제범주 \(y_i\)	\(\pi_{i1}\)	\(\pi_{i2}\)	\(\pi_{i3}\)	추정범주 \(\hat{y}_i\)
25	5	3	0.0026	0.0107	0.9867	3
25	10	3	0.0116	0.0452	0.9432	3
25	20	2	0.1905	0.3567	0.4528	3
25	30	1	0.8252	0.1352	0.0396	1
32	5	3	0.0124	0.0481	0.9395	3
32	10	3	0.0531	0.1706	0.7763	3
32	20	2	0.5296	0.3230	0.1474	1
32	30	1	0.9576	0.0338	0.0085	1
42	5	1	0.1049	0.2709	0.6241	3
42	10	3	0.3443	0.3852	0.2705	2
42	20	1	0.9133	0.0685	0.0181	1
42	30	1	0.9953	0.0038	0.0009	1

6.4.3 인근범주 로짓모형

인근범주 로짓모형은 아래와 같이 인접한 두 범주의 확률 비율에 대한 회귀모형이다.

\[\begin{equation} \log \left( \frac{\pi_{ij}}{\pi_{i,j+1}} \right) = \beta_{0,j} + \boldsymbol\beta^\top \mathbf{x}_i, \, j = 1, \cdots, J - 1 \tag{6.11} \end{equation}\]

따라서, \(\pi_{ij}\)간에 다음과 같은 관계식이 성립한다.

\[\begin{equation*} \begin{split} \pi_{ij} &= \exp (\beta_{0,j} + \boldsymbol\beta^\top \mathbf{x}_i) \pi_{i,j+1}\\ &= \pi_{iJ} \exp \left(\sum_{k = j}^{J - 1} \beta_{0,k} + (J - j) \boldsymbol\beta^\top \mathbf{x}_i\right), \, j = 1, \cdots, J-1\\ \sum_{j = 1}^{J} \pi_{ij} &= 1 \end{split} \end{equation*}\]

이를 정리하면

\[\begin{equation} \begin{split} \pi_{ij} &= \frac{\exp \left( \sum_{l = j}^{J - 1} \beta_{0,l} + (J - j) \boldsymbol\beta^\top \mathbf{x}_i \right)}{1 + \sum_{k = 1}^{J - 1} \exp \left( \sum_{l = k}^{J - 1} \beta_{0,l} + (J - k) \boldsymbol\beta^\top \mathbf{x}_i \right)}, \, j = 1, \cdots, J - 1\\ \pi_{iJ} &= \frac{1}{1 + \sum_{k = 1}^{J - 1} \exp \left( \sum_{l = k}^{J - 1} \beta_{0,l} + (J - k) \boldsymbol\beta^\top \mathbf{x}_i \right)} \end{split} \tag{6.12} \end{equation}\]

와 같다. 이는 6.3.2절에서 살펴보았던 명목형 로지스틱 회귀모형에 비해 다소 복잡하지만 비슷한 형태이며, 역시 최우추정법을 이용하여 모형을 추정할 수 있다.

R에서는 VGAM 패키지의 vglm 함수를 이용할 때 파라미터 family의 값을 VGAM 패키지의 acat 함수를 설정함으로써 인근범주 로짓모형을 추정할 수 있다. 이 때 acat 함수의 parallel 파라미터값을 TRUE로 설정함으로써 식 (6.11)에서와 같이 비례 승산 모형(proportional odds model)을 정의한다.

vglm_fit <- VGAM::vglm(
  y ~ N + L,
  data = train_df,
  family = VGAM::acat(reverse = TRUE, parallel = TRUE)
  )

print(coef(vglm_fit))

## (Intercept):1 (Intercept):2             N             L 
##    -9.0658976    -8.9018134     0.1725867     0.2082167

추정된 모형을 위 식 (6.12)에 대입하면 각 객체 \(i\)가 각 범주 \(j\)에 속할 확률을 추정할 수 있다. VGAM 패키지의 predictvglm 함수가 해당 계산을 수행한다.

predict_df <- VGAM::predictvglm(vglm_fit, train_df, "response") %>% 
  as_data_frame() %>%
  `colnames<-`(c("p1", "p2", "p3")) %>%
  mutate(pred_class = ordered(apply(., 1, function(x) which.max(x))))

bind_cols(train_df, predict_df) %>%
  knitr::kable(
    booktabs = TRUE,
    align = c('r', 'r', 'r', 'r', 'r', 'r', 'r'),
    col.names = c('잡음(N)', '손실(L)', '실제범주 $y_i$', 
                  '$\\pi_{i1}$', '$\\pi_{i2}$', '$\\pi_{i3}$', '추정범주 $\\hat{y}_i$'),
    caption = '위 통신 만족도 데이터에 대한 인근범주 로짓모형의 추정범주',
#    caption = 'Table \\@ref(tab:ordinal-logistic-reg-train-data)에 대한 인근범주 로짓모형의 추정범주',
    digits = 4
  )

Table 6.8: 위 통신 만족도 데이터에 대한 인근범주 로짓모형의 추정범주

잡음(N)	손실(L)	실제범주 \(y_i\)	\(\pi_{i1}\)	\(\pi_{i2}\)	\(\pi_{i3}\)	추정범주 \(\hat{y}_i\)
25	5	3	0.0007	0.0280	0.9713	3
25	10	3	0.0052	0.0751	0.9197	3
25	20	2	0.1804	0.3244	0.4952	3
25	30	1	0.7894	0.1770	0.0337	1
32	5	3	0.0072	0.0874	0.9054	3
32	10	3	0.0474	0.2045	0.7480	3
32	20	2	0.5611	0.3015	0.1375	1
32	30	1	0.9339	0.0626	0.0036	1
42	5	1	0.1393	0.3026	0.5581	3
42	10	3	0.4411	0.3385	0.2204	1
42	20	1	0.9063	0.0867	0.0070	1
42	30	1	0.9881	0.0118	0.0001	1

비례 승산 모형(proportional odds model)이 아닌 인근범주 로짓모형은 아래와 같다.

\[\begin{equation} \log \left( \frac{\pi_{ij}}{\pi_{i,j+1}} \right) = \beta_{0,j} + \boldsymbol\beta_j^\top \mathbf{x}_i, \, j = 1, \cdots, J - 1 \tag{6.13} \end{equation}\]

해당 모형은 acat 함수의 parallel 파라미터 값을 FALSE로 설정함으로써 추정할 수 있다.

vglm_fit <- VGAM::vglm(
  y ~ N + L,
  data = train_df,
  family = VGAM::acat(reverse = TRUE, parallel = FALSE)
  )

print(coef(vglm_fit))

## (Intercept):1 (Intercept):2           N:1           N:2 
##  -12.94227208   -6.10597713    0.31431599    0.04643039 
##           L:1           L:2 
##    0.17500408    0.29578067

predict_df <- VGAM::predictvglm(vglm_fit, train_df, "response") %>% 
  as_data_frame() %>%
  `colnames<-`(c("p1", "p2", "p3")) %>%
  mutate(pred_class = ordered(apply(., 1, function(x) which.max(x))))

bind_cols(train_df, predict_df) %>%
  knitr::kable(
    booktabs = TRUE,
    align = c('r', 'r', 'r', 'r', 'r', 'r', 'r'),
    col.names = c('잡음(N)', '손실(L)', '실제범주 $y_i$', 
                  '$\\pi_{i1}$', '$\\pi_{i2}$', '$\\pi_{i3}$', '추정범주 $\\hat{y}_i$'),
    caption = '위 통신 만족도 데이터에 대한 인근범주 로짓모형의 추정범주 (비례 승산 모형이 아닌 경우)',
#    caption = 'Table \\@ref(tab:ordinal-logistic-reg-train-data)에 대한 인근범주 로짓모형의 추정범주 (비례 승산 모형이 아닌 경우)',
    digits = 4
  )

Table 6.9: 위 통신 만족도 데이터에 대한 인근범주 로짓모형의 추정범주 (비례 승산 모형이 아닌 경우)

잡음(N)	손실(L)	실제범주 \(y_i\)	\(\pi_{i1}\)	\(\pi_{i2}\)	\(\pi_{i3}\)	추정범주 \(\hat{y}_i\)
25	5	3	0.0004	0.0303	0.9693	3
25	10	3	0.0043	0.1200	0.8757	3
25	20	2	0.1295	0.6313	0.2392	2
25	30	1	0.5365	0.4546	0.0089	1
32	5	3	0.0055	0.0412	0.9533	3
32	10	3	0.0488	0.1517	0.7995	3
32	20	2	0.5924	0.3200	0.0876	1
32	30	1	0.9131	0.0857	0.0012	1
42	5	1	0.1667	0.0536	0.7797	3
42	10	3	0.6335	0.0850	0.2815	1
42	20	1	0.9734	0.0227	0.0039	1
42	30	1	0.9959	0.0040	0.0000	1

References

Cavanaugh, JR, RW Hatch, and JL Sullivan. 1976. “Models for the Subjective Effects of Loss, Noise, and Talker Echo on Telephone Connections.” Bell System Technical Journal 55 (9): 1319–71.

Czepiel, Scott A. 2002. “Maximum Likelihood Estimation of Logistic Regression Models: Theory and Implementation.” Available at Czep. Net/Stat/Mlelr. Pdf.

Pratt, John W. 1981. “Concavity of the Log Likelihood.” Journal of the American Statistical Association 76 (373): 103–6.

전치혁. 2012. 데이터마이닝 기법과 응용. 한나래출판사.