Non-linear models

앞서 언급했다시피 이항변수일 경우에도 우리는 OLS를 통해 예측변수와 종속변수 간 관계를 추정할 수 있습니다. 이를 선형확률모형(Linear Probability Models, LPM)이라고 하는데, 이 경우 추정된 선이 0-1을 넘어서는 외삽의 문제(extrapolate problem)가 제기될 수 있습니다.

한편, Logit 모델은 종속변수가 이항변수일 때 사용하는 비선형 모델입니다.

  • 즉, 우리가 관심을 가지고 이해하고자 하는, 혹은 모형화하고자 하는 변수가 두 개의 값을 취한다는 것을 의미합니다.

  • 성공일 경우에는 1, 실패일 경우에는 0의 값을 갖습니다.

    • 예를 들어, 어떤 예측변수가 개인의 결혼여부에 영향을 미치는지 알고싶어한다고 해보겠습니다. 이때, 기혼은 1로, 미혼은 0인 두 개의 값만을 가지게 될 것입니다.

    • 그러나 성패를 어떻게 개념정의할 것인지는 전적으로 연구자에게 달려있습니다.

      • 과부의 경우 결혼 상태로 볼 것인가 혹은 미혼 상태로 볼까요 혹은 별거 중인 사람은 어떻게 판단해야 할까요?

      • 6년 동안 약혼하고 함께 살았으면 결혼한 것으로 간주해야할까요? 아니면 법적으로 결혼한 것이 아니니 미혼으로 처리해야 할까요?

    • 우리가 단순하게 법적 결혼에 대해 이해하고자 하는 것이 아니라 결혼을 어떤 개인이 장기에 걸쳐 맺기를 원하는 관계라는 측면에서 이해하고자 한다면, 1과 0으로 소급하는 자료 유형은 연구자의 판단을 필요로 합니다.

  • 개념화의 문제는 전적으로 연구자에게 달려있지만 모델은 오직 0과 1의 값을 갖는 더미변수의 이항변수의 형태를 가질 것이 선호됩니다.

또, Logit 모델의 경우는 비선형성을 가정하기 때문에 OLS와 마찬가지로 편미분을 통해 얻어지는 계수값을 확률로 직접 해석하는 데 무리가 있습니다. 따라서 Logit 분석을 통해 얻은 계수를 해석하는 데에는 두 가지 접근법이 있습니다.

  1. MEM (Marginal Effect at the Mean)

  • 변수들을 평균에 고정시키고 보고자 하는 에측변수의 한 단위 증가가 종속변수에 가져오는 변화를 보는 방법입니다.

  • \(x\)의 평균을 기준으로 단위를 변화시켜가며 그에 따른 \(y\)의 차이를 통해 한계효과를 확인하는 방법입니다.

  1. AME (Average Marginal Effects)

  • \(x\)의 각 관측된 값들에 따른 \(y\)의 값을 구하고, 각 \(x\)에 일정한 단위를 증가시킵니다.

  • 그러면 우리는 각 \(x_i\)와 \(x_{i+1}\)에 따른 \(y\) 차이를 i개 만큼 갖게 될 것입니다.

  • 그 \(i\)개의 차이들의 평균을 구하는 것입니다.

그럼 승산비(odds ratio)란 무엇일까요? \(\text{Odds} = \frac{P}{1-P}\).

  • 경기에서 이길 확률이 1/5, 질 확률이 4/5라면 게임에서 이길 승산비는 1/4가 되며, 계산된 값을 바탕으로 우리는 5번 중에 4번 질 동안 1번 이긴다라고 해석할 수 있습니다.

    • 승산비가 1보다 크면 예측 변수가 증가함에 따라 사건 발생 확률이 증가한다는 것을 나타냅니다.

  • 승산비가 1보다 작으면 예측 변수가 증가함에 따라 사건 발생 확률이 감소한다는 것을 나타냅니다.

일단 R을 통해 예제데이터를 살펴보며 익숙해져 보도록 하겠습니다.

필요한 패키지를 로드

library(curl)
library(MASS)
library(BBmisc)
library(here)
library(ezpickr)
library(broom)
library(tidyverse) 
library(showtext)
library(extrafont)
loadfonts()
font_add("NanumGothic", "NanumGothic.ttf")
theme_set(
  theme_bw() +
    theme(text = element_text(family = "NanumGothic"))
)
Sys.setlocale("LC_CTYPE","korean")
## [1] "Korean_Korea.949"
## 데이터를 불러옵니다.
df <- pick("heart_data.dta")
glimpse(df)
## Rows: 5,388
## Columns: 12
## $ death    <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0...
## $ anterior <dbl+lbl>  1,  0,  0,  1,  1,  1,  0, NA, NA,  0,  1, NA,  0,  0...
## $ hcabg    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ kk1      <dbl> 1, 0, 1, 0, 0, NA, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, NA,...
## $ kk2      <dbl> 0, 1, 0, 1, 1, NA, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, NA,...
## $ kk3      <dbl> 0, 0, 0, 0, 0, NA, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA,...
## $ kk4      <dbl> 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA,...
## $ age1     <dbl> 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0...
## $ age2     <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0...
## $ age3     <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1...
## $ age4     <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ center   <dbl> 1255, 1255, 1255, 1255, 1255, 1255, 1255, 1255, 1255, 1255...

불러온 heart_data.dta는 STATA 파일입니다.

R에서 Logit 모델 추정

먼저 Logit 모델을 추정해보겠습니다. OLS와 비슷하지만 Logit 모델은 최소자승법(ordinary least square, OLS)이 아니라 일반화 선형 모형(generalized linear model, GLM)을 사용합니다.

model <- glm(
  death ~ anterior + hcabg + kk2 + kk3 + kk4 + 
    age2 + age3 + age4, 
  data = df,
  family = binomial(link = "logit"))

여기서 절편(intercept)은 모든 예측변수들이 0일 때의 예측된 종속변수의 값입니다. 우리가 관심이 있는 것은 예측변수와 종속변수 간의 관계이기 때문에 절편은 의미가 없습니다. 따라서 일반적으로 해석하지 않습니다.

우선 모형의 결과를 살펴보겠습니다.

summary(model)
## 
## Call:
## glm(formula = death ~ anterior + hcabg + kk2 + kk3 + kk4 + age2 + 
##     age3 + age4, family = binomial(link = "logit"), data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3670  -0.3276  -0.1986  -0.1443   3.1807  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.0521     0.2586 -19.536  < 2e-16 ***
## anterior      0.6426     0.1676   3.835 0.000126 ***
## hcabg         0.7444     0.3530   2.109 0.034929 *  
## kk2           0.8117     0.1805   4.497 6.90e-06 ***
## kk3           0.7757     0.2691   2.883 0.003939 ** 
## kk4           2.6597     0.3560   7.471 7.96e-14 ***
## age2          0.4930     0.3102   1.589 0.111961    
## age3          1.5112     0.2662   5.676 1.38e-08 ***
## age4          2.1853     0.2718   8.039 9.06e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1484.6  on 4482  degrees of freedom
## Residual deviance: 1273.3  on 4474  degrees of freedom
##   (905 observations deleted due to missingness)
## AIC: 1291.3
## 
## Number of Fisher Scoring iterations: 7

summary()를 통해 나타나는 분석결과는 약간 보기가 난잡하니 분석이 용이하도록 그 결과를 tibble() 형태로 바꾸어주도록 하겠습니다.

tidy_model <- model %>% tidy(); tidy_model

이번에는 신뢰구간을 구해보도록 하겠습니다.

tidy_model <- tidy_model %>%  # 이미 만들어놓은 Logit 결과 tibble에
  bind_cols( # 변수를 하나 추가하는데,
    model %>% broom::confint_tidy() %>% # 그 변수는 broom 패키지의 신뢰구간을 계산해주는
      ## confint_tidy() 함수를 이용
      select(ll = conf.low, ul = conf.high) # 그렇게 만들어질 변수 이름을 ll, ul로 바꾸어 저장
  ); tidy_model
# 변수들의 이름에 들어가 있는 '.'를 '_'로 변경
tidy_model <- tidy_model %>%
  rename_at(vars(contains('.')), function(x) str_replace(x, '\\.', '_')); tidy_model
## regular expression에서 그냥 .은 말그대로 '아무거나'(anything)을 의미
## 정규식에서의 .과 우리가 지정하고 싶은 '.'을 구분해주기 위해서 \\.라고 지정

모델로부터 예측값과 예측확률을 구하기

death의 확률은 위에서 구한 tidy_model의 estimate로 재구성한 Logit 회귀식에 각 관측치 값을 투입함으로써 계산할 수 있습니다. 혹은 predict() 함수나 broom 패키지의 augment() 함수를 이용해서도 계산할 수 있습니다.

data_augmented <- broom::augment(model); glimpse(data_augmented)
## Rows: 4,483
## Columns: 16
## $ .rownames  <chr> "1", "2", "3", "4", "5", "7", "10", "11", "13", "14", "1...
## $ death      <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,...
## $ anterior   <dbl+lbl> 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ hcabg      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ kk2        <dbl> 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,...
## $ kk3        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ kk4        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ age2       <dbl> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,...
## $ age3       <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0,...
## $ age4       <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ .fitted    <dbl> -4.409516, -2.055082, -3.540901, -2.086646, -3.104792, -...
## $ .resid     <dbl> -0.1554845, 2.0859534, -0.2390607, 2.0993481, -0.2961685...
## $ .std.resid <dbl> -0.1555426, 2.0907571, -0.2391672, 2.1027987, -0.2965342...
## $ .hat       <dbl> 0.0007472251, 0.0045899510, 0.0008899623, 0.0032792762, ...
## $ .sigma     <dbl> 0.5335234, 0.5326118, 0.5335165, 0.5326012, 0.5335100, 0...
## $ .cooksd    <dbl> 1.011183e-06, 4.018575e-03, 2.871499e-06, 2.955340e-03, ...

data_augmented는 우리가 Logit 분석으로 구한 계수, estimate를 이용해 회귀식을 재구성한 뒤, 개별 관측치들을 그 식에 하나하나 투입하여 death에 대한 예측값을 구한 것입니다. 여기서는 연습이니까 예측값을 구하는 get_fitted_values()라는 함수를 직관적으로 구성해 보겠습니다.

get_fitted_values <- function() {
  
  ## 모델에서 사용할 변수들만 선별
  df <- df %>% select(anterior, hcabg, kk2:kk4, age2:age4)
  ## Logit 분석에서 구한 estimates를 저장해주되, 절편값(intercept)는 제외
  estimates <- tidy_model$estimate[2:length(tidy_model$estimate)]
  ## 결과적으로 우리는 하나의 행에 선형예측변수들을 담을 수 있습니다.
  ## 그리고 이 과정은 개별 관측치별로 반복하는 루프문을 구성
  get_fitted_value <- function(x) {
    row_vector <- df %>% slice(x) %>% as.numeric()
    result <- row_vector * estimates
    result <- sum(c(result, tidy_model$estimate[1])) 
    ## 만약 NA가 있으면, 그냥 NA로 리턴
    return(result)
  }
  ## 이 루프문이 각 행마다 반복될 수 있도록 functional을 구성.
  ## 여기서는 double 유형이기 때문에 map_dbl()이라는 functional을 사용
  model_fitted_values <- 1:nrow(df) %>% map_dbl(get_fitted_value)
  # 결과를 리턴
  return(model_fitted_values)
}
fitted_values <- get_fitted_values()
fitted_values[1:10]
##  [1] -4.409516 -2.055082 -3.540901 -2.086646 -3.104792        NA -5.052071
##  [8]        NA        NA -4.559047
data_augmented$.fitted[1:10]
##  [1] -4.409516 -2.055082 -3.540901 -2.086646 -3.104792 -5.052071 -4.559047
##  [8] -4.409516 -5.052071 -3.540901

get_fitted_values()는 augment() 결과와 약간 차이가 있을 수도 있습니다. 왜냐하면 augment()는 자동으로 모델에 포함되는 사례들만 계산해주기 때문입니다 (표본 추정). 우리가 직접 만든 함수는 NA가 있을 때, NA를 반환하기는 하지만 Logit 모델에 포함되는 값들을 선별해주는 작업은 제외되어 있습니다.

new_data <- df %>% 
  dplyr::select(anterior, hcabg, kk2:kk4, age2:age4) %>% drop_na()

predict(model, newdata = new_data)[1:10]
##         1         2         3         4         5         6         7         8 
## -4.409516 -2.055082 -3.540901 -2.086646 -3.104792 -5.052071 -4.559047 -4.409516 
##         9        10 
## -5.052071 -3.540901

하지만 예측값은 그 자체로 아무것도 전달하지 않습니다. 따라서 우리는 예측값들로부터 death에 대한 예측확률(predicted probability)를 구하고자 합니다. Logit 계수값으로부터 확률을 얻는 공식은 다음과 같습니다: \(\text{Pr(logit }\theta)= \frac{1}{1+ e^{-x}}\). 이번에도 연습은 연습이니까 get_probability()라는 예측확률을 구하는 함수를 get_fitted_values()의 결과를 활용하여 만들어 보겠습니다.

get_probablity <- function(x) {
  
  if (is.na(x)) { ## 결측치인지 확인
    return(NA_real_) # x가 결측치면, real 유형의 NA를 리턴!
  }
## Logit 계수값을 지수화하는 함수
  result <- 1 / (1 + exp(-x))
  return(result)
}
get_predicted_probabilities <- function() {
  result <- map_dbl(fitted_values, get_probablity) #functional을 이용
  return(result) # fitted_values로부터 예측확률을 계산
}
our_predictions <- get_predicted_probabilities()
new_data <- df %>% select(anterior, hcabg, kk2:kk4, age2:age4)
predictions_r <- predict(model, 
                         newdata = new_data, 
                         type = 'response') %>% unname()
## 물론 predict()에 type 옵션을 이용하는 것이 훨씬 효율적
our_predictions[1:10]
##  [1] 0.012014954 0.113539925 0.028170609 0.110401541 0.042910041          NA
##  [7] 0.006355425          NA          NA 0.010363512
predictions_r[1:10]
##  [1] 0.012014954 0.113539925 0.028170609 0.110401541 0.042910041          NA
##  [7] 0.006355425          NA          NA 0.010363512
## 두 결과가 같다는 것을 확인

Logit 모델의 적합도 (Fitness of the Logit model)

모델은 종속변수와 예측변수 간의 관계를 일반적으로 설명하기 위하여 사용된다고 할 수 있습니다.

  • 따라서 모델의 본연의 목적은 개별 관측치를 완벽하게 예측하는 모델을 구축하는 것은 아닙니다(과적합의 문제, overfitting).

  • 만약 데이터가 모델과 지나치게 잘 맞아떨어진다면, 우리는 표본을 넘어서 모집단에 속하는 사람들의 death의 확률을 그 모델로 잘 예측하지 못할 수 있습니다.

  • 즉, “표본은 완벽하게 예측하는 모델”일지는 몰라도 그것이 “모집단을 추론하는 데 도움을 주는 모델”은 아닐 수 있다는 것입니다.

  • 만약 모델이 과적합이라면 표본을 통해서는 death를 야기할 수 있는 요인들을 이해하는 데 도움을 줄 수 있지만, 정작 모집단에 일반화할 수는 없는, 그런 문제가 발생할 수 있습니다.

모델의 적합도를 측정하는 한 가지 측정지표는 바로 혼동행렬(confusion matrix)입니다. 혼동행렬은 얼마나 예측이 잘 수행되는지, 즉 얼마나 실제 원데이터의 값을 모델로 알 수 있는지를 보여줍니다.

library(regclass)
model_confusion_matrix <- confusion_matrix(model)
model_confusion_matrix
##          Predicted 0 Predicted 1 Total
## Actual 0        4300           7  4307
## Actual 1         172           4   176
## Total           4472          11  4483
## 직접 매뉴얼대로 기본 함수를 이용해서 confusion matrix를 구할 수도 있지만, 
## 여기서는 패키지를 이용해보겠습니다. regclass 패키지
## Logit 모델을 추정하여 저장한 augmented 데이터에 예측확률을 추가
data_augmented <- data_augmented %>%
  mutate( # 예측확률과 예측 승산(odds)를 계산하는데, functional, map_dbl()을 이용해
    ## 이전에 만들어 두었던 get_probability() 함수를 반복적용
    .predicted_probabilities = map_dbl(.fitted, get_probablity),
    .predicted_odds = .predicted_probabilities / 
      (1 - .predicted_probabilities)
  )

data_augmented %>% glimpse()
## Rows: 4,483
## Columns: 18
## $ .rownames                <chr> "1", "2", "3", "4", "5", "7", "10", "11", ...
## $ death                    <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ anterior                 <dbl+lbl> 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0,...
## $ hcabg                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ kk2                      <dbl> 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, ...
## $ kk3                      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ kk4                      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ age2                     <dbl> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, ...
## $ age3                     <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, ...
## $ age4                     <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ .fitted                  <dbl> -4.409516, -2.055082, -3.540901, -2.086646...
## $ .resid                   <dbl> -0.1554845, 2.0859534, -0.2390607, 2.09934...
## $ .std.resid               <dbl> -0.1555426, 2.0907571, -0.2391672, 2.10279...
## $ .hat                     <dbl> 0.0007472251, 0.0045899510, 0.0008899623, ...
## $ .sigma                   <dbl> 0.5335234, 0.5326118, 0.5335165, 0.5326012...
## $ .cooksd                  <dbl> 1.011183e-06, 4.018575e-03, 2.871499e-06, ...
## $ .predicted_probabilities <dbl> 0.012014954, 0.113539925, 0.028170609, 0.1...
## $ .predicted_odds          <dbl> 0.012161069, 0.128082390, 0.028987196, 0.1...
## 이제 예측 대 실제(predicted_vs_actual)라는 이름의 새 변수를 만들어 주자.

data_augmented <- data_augmented %>%
  mutate(
    predicted_vs_actual = case_when(
      death == 0 & .predicted_probabilities < 0.5 ~ 'raw0_pred0',
      death == 0 & .predicted_probabilities >= 0.5 ~ 'raw0_pred1',
      death == 1 & .predicted_probabilities < 0.5 ~ 'raw1_pred0',
      death == 1 & .predicted_probabilities >= 0.5 ~ 'raw1_pred1',
      T ~ NA_character_
    )
  )

case_when()을 이용하여 실제 death의 값이 0일 때와 우리가 구한 예측확률이 0.5 미만일 때와 같이 실제 death의 값(2)과 예측확률을 0.5 기준으로 구간화하여(2) \(2\times 2\) 교차표를 만들어주는 것입니다.

## 이렇게 만든 예측 대 실제 변수를 요인형(factor)로 변환
data_augmented <- data_augmented %>%
  mutate(
    predicted_vs_actual = parse_factor(
      predicted_vs_actual, levels = c(
        'raw0_pred0', 'raw0_pred1', 'raw1_pred0', 'raw1_pred1'
      )
    )
  )
glimpse(data_augmented)
## Rows: 4,483
## Columns: 19
## $ .rownames                <chr> "1", "2", "3", "4", "5", "7", "10", "11", ...
## $ death                    <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ anterior                 <dbl+lbl> 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0,...
## $ hcabg                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ kk2                      <dbl> 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, ...
## $ kk3                      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ kk4                      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ age2                     <dbl> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, ...
## $ age3                     <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, ...
## $ age4                     <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ .fitted                  <dbl> -4.409516, -2.055082, -3.540901, -2.086646...
## $ .resid                   <dbl> -0.1554845, 2.0859534, -0.2390607, 2.09934...
## $ .std.resid               <dbl> -0.1555426, 2.0907571, -0.2391672, 2.10279...
## $ .hat                     <dbl> 0.0007472251, 0.0045899510, 0.0008899623, ...
## $ .sigma                   <dbl> 0.5335234, 0.5326118, 0.5335165, 0.5326012...
## $ .cooksd                  <dbl> 1.011183e-06, 4.018575e-03, 2.871499e-06, ...
## $ .predicted_probabilities <dbl> 0.012014954, 0.113539925, 0.028170609, 0.1...
## $ .predicted_odds          <dbl> 0.012161069, 0.128082390, 0.028987196, 0.1...
## $ predicted_vs_actual      <fct> raw0_pred0, raw1_pred0, raw0_pred0, raw1_p...
table(data_augmented$predicted_vs_actual)
## 
## raw0_pred0 raw0_pred1 raw1_pred0 raw1_pred1 
##       4300          7        172          4

이번에는 예측확률이 얼마나 잘 맞아 떨어졌는지를 계산해보겠습니다. 먼저 예측 대 실제 변수가 가지는 결측치를 제외한 모든 관측치의 수를 total_cases라는 객체에 저장해줍니다. 그리고 4개의 카테고리가 갖는 빈도값을 하나의 벡터로 저장해줍니다. 그리고 예측과 실제가 맞아떨어진 경우(0-0, 1-1)이 전체 실제 사례 중 몇 퍼센트를 차지하는지 계산하는 것입니다.

total_cases <- table(data_augmented$predicted_vs_actual) %>%
  as.numeric() %>% sum()
table_vector <- table(data_augmented$predicted_vs_actual) %>% as.numeric()
percent_correct <- (table_vector[1] + table_vector[4]) / total_cases
percent_correct   # 96% 의 사례를 예측해낸 것을 확인
## [1] 0.9600714

Logit 모형의 적합도를 확인할 수 있는 또 다른 측정지표는 반응자 작용특성(Reciever Operating Characteristic, ROC) 곡선입니다. 적중확률(Y축, True positive rate, sensityvity) 대 오경보확률(X축, False positive rate, 1-Specificity)를 그래프로 나타내는 것입니다. 큰 값을 가질수록 적합도가 높은 것으로 판단할 수 있습니다.

## pROC라는 패키지
library(pROC)
data_for_roc <- data_augmented # logit 모형에 사용될 데이터를 augmented
model_for_roc <- glm(
  death ~ anterior + hcabg + kk2 + kk3 + kk4 + age2 + age3 + age4, 
  data = data_for_roc, family = binomial(link = 'logit')
)
predictions_for_roc <- predict(model_for_roc, type = "response")
## 개별 관측치에 대한 예측확률을 계산

data_for_roc$prob <- predictions_for_roc # 이렇게 구한 예측확률을 prob이란 변수로
roc_curve <- roc(death ~ prob, data = data_for_roc) 
## 그리고 이항종속변수와 이 예측확률 간 관계를 계산한 반응자 작용특성, 
## ROC를 roc_curve라고 하는 데이터에 저장
plot(roc_curve)

## 직선에서 멀리 벗어날수록 모델 적합도가 높다고 할 수 있습니다.

혼동행렬과 ROC 말고 실제값과 예측값을 비교하여 보여줄 수 있는 적합도에 관한 통계치는 바로 Hosmer-Lemeshow의 적합도 (Goodness of Fit, GOF) 통계치입니다.

library(ResourceSelection) # 새로운 패키지를 로드
data_for_lmgof <- data_augmented
model_for_lmgof <- glm(
  death ~ anterior + hcabg + kk2 + kk3 + kk4 + age2 + age3 + age4, 
  data = data_for_lmgof, family = binomial(link = 'logit')
)
hoslem.test(data_for_lmgof$death, fitted(model_for_lmgof), g = 10)
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  data_for_lmgof$death, fitted(model_for_lmgof)
## X-squared = 9.4776, df = 8, p-value = 0.3036
## GOF 통계치 값이 클수록 작은 p-values를 가지게 되는데, 이는 적합도가 낮다는 것
## 높은 p-value를 가지는 경우, 적합도가 높다는 것을 의미

모수 추정의 해석

보통 OLS의 경우에는 계수 그 자체를 가지고 부호, 유의도, 그리고 크기를 다른 계수들과 비교하여 해석할 수 있습니다.

  • 하지만 Logit 모델의 계수 그 자체는 바로 해석하기 까다롭습니다.

  • 왜냐하면 원래는 비선형 관계이기 때문에 \(x\)가 0에서 1로 증가했을 때와 1에서 2로 증가했을 때의 변화분이 같지 않으나, 그것을 로그변환을 통해 강제로 선형화시켜준 결과이기 때문입니다.

    • 따라서 Logit의 계수값은 그 자체로는 확률적으로 예측변수와 종속변수 간의 관계를 설명해주지 못합니다.

    • 우리가 앞서 복잡한 함수들을 이용하여 예측확률을 구해주었던 이유입니다.

  • 그러나 예측확률은 그래프 등을 통하여 예측변수와 종속변수 간의 관계 추이를 보여줄 수는 있지만 그 둘의 관계를 요약하여 보여줄 수는 없습니다.

여기서 한계효과(marginal effect)가 등장합니다.

tidy_model # tidy하게 정리한 Logit 결과를 다시 불러와 보겠습니다.

선형회귀분석의 계수값들이 모두 같은 방식으로 해석될 수 있다는 것은 잘 알고 있는 사실입니다.

  • 선형회귀분석의 계수값들은 예측변수의 한 단위 변화가 종속변수에 얼마만큼의 변화와 연관되는가를 보여줍니다.

  • Logit 모델에서 승산비(Odds ratios)는 선형회귀분석의 계수값과 같은 방식으로 해석할 수 있습니다.

  • 승산비를 구하기 위해서는 로그화된 Logit 계수값들을 지수화(exponentiate)해야 합니다.

exp(tidy_model$estimate)
## [1]  0.006396075  1.901332984  2.105274605  2.251732160  2.172105099
## [6] 14.291372300  1.637260198  4.532028652  8.893222052
## 계수값 뿐 아니라 신뢰구간도 지수화될 수 있습니다.
exp(tidy_model$ll)
## [1] 0.003731944 1.373610083 0.994274803 1.574745487 1.251910232 6.991113190
## [7] 0.894460373 2.744267109 5.319540393
exp(tidy_model$ul)
## [1]  0.01032747  2.65232006  4.02161372  3.19969131  3.61189219 28.42266646
## [7]  3.04323642  7.83821034 15.53112158
## mfx 패키지를 이용해서 승산비를 구할 수 있습니다.
library(mfx)
logit_or_model <- logitor(model, data = data_augmented)
logit_or <- logit_or_model$oddsratio %>% as.data.frame() %>% as_tibble() %>%
  tibble::rowid_to_column() %>% select(
    odds_ratio = OddsRatio, se = `Std. Err.`, z, p_value = `P>|z|`
  ) # 복잡해보이지만 간단합니다. 
## Logit 결과값의 계수를 승산비로 바꾸어주어 
## logit_or_model에 저장하고, 
## 이를 tibble로 바꾸어준 다음에, 행의 라인을 변수화 시켜
## tibble을 정리해준 이후에 필요한 변수들을 이름을 바꾸어 
## select로 뽑아 저장한 것입니다.
logit_or

mfx 패키지로 손쉽게 Logit 계수값을 승산비로 바꾸어줄 수 있지만, 이 패키지는 주로 의미있게 해석되지는 않지만 모형에 포함할 필요는 있는 ’절편값’을 날려버립니다. 따라서 절편값은 직접 계산해서 변수로 추가해주어야 합니다. 혹은 아예 손수 처음부터 tidy_model에서 logit_or과 같은 tidy_model_or을 구할 수도 있습니다.

tidy_model_or <- tidy_model %>% 
  mutate(
    odds_ratio = exp(estimate),  # 승산비/기울기
    var.diag = diag(vcov(model)),  # 각 계수값의 분산
    or.se = sqrt(odds_ratio^2 * var.diag),  # 계수값의 분산으로 조정된 승산비의 표준오차
    statistic = qnorm(p_value / 2), # z값
    ll = exp(ll), ul = exp(ul) # 신뢰구간을 지수화
  ) %>%
  select(term, odds_ratio, se = or.se, statistic, p_value, ll, ul)
tidy_model_or

이번에는 Logit 계수값을 표준화해보겠습니다.

  • 표준화를 시키면 예측변수들 간의 계수 크기의 비교가 가능해집니다.

    • 별로 추천하고 싶은 방법은 아닙니다. 어차피 사회과학에서 특정 변수의 효과가 다른 변수보다 종속변수에 미치는 효과가 크다라고 주장하는 것은 확률적으로 표본에 따라 반증될 가능성도 높고, 통계적으로 그러하다고 해서 실제로 그러하리라고 기대할만한 메커니즘을 찾기도 힘들기 때문입니다.

  • 여하튼 표준화하는 방법은 다음의 공식을 따릅니다: \(\text{standardaized} = \frac{(x - \mu)}{\sigma}\).

library(reghelper) # reghelper라는 패키지를 이용
reghelper::beta(model)
## 
## Call:
## glm(formula = "death ~ anterior.z + hcabg.z + kk2.z + kk3.z + kk4.z + age2.z + age3.z + age4.z", 
##     family = binomial(link = "logit"), data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3670  -0.3276  -0.1986  -0.1443   3.1807  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.72682    0.11354 -32.824  < 2e-16 ***
## anterior.z   0.31980    0.08339   3.835 0.000126 ***
## hcabg.z      0.12950    0.06140   2.109 0.034929 *  
## kk2.z        0.32341    0.07192   4.497 6.90e-06 ***
## kk3.z        0.17150    0.05949   2.883 0.003939 ** 
## kk4.z        0.27376    0.03664   7.471 7.96e-14 ***
## age2.z       0.21661    0.13628   1.589 0.111961    
## age3.z       0.66152    0.11654   5.676 1.38e-08 ***
## age4.z       0.71194    0.08856   8.039 9.06e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1484.6  on 4482  degrees of freedom
## Residual deviance: 1273.3  on 4474  degrees of freedom
## AIC: 1291.3
## 
## Number of Fisher Scoring iterations: 7
## reghelper의 beta()라는 함수는 모든 계수값을 표준화

이 결과는 나이가 더 많은 집단일수록 종속변수를 설명하는 힘이 큰 예측변수라는 것을 확인할 수 있습니다.

한계효과(marginal effect)

한계효과는 예측변수의 한 단위 변화가 종속변수에 가져오는 효과를 의미합니다. 그러나 선형회귀분석의 계수값처럼 비선형회귀분석에서의 승산비는 일정하게 나타나지가 않기 때문에 그 변화를 어떻게 보여주는가가 중요합니다.

library(margins)
## margins 패키지의 함수로 한계효과를 구합니다.
marginal_effects_from_model <- margins(model) 
## 각 변수의 평균 한계효과와 그 표준오차, 유의수준, 신뢰구간 등을 확인할 수 있습니다.
summary(marginal_effects_from_model) 
marginal_effects_df <- summary(marginal_effects_from_model) %>%
  as.data.frame() %>% as_tibble()
## 이렇게 구한 한계효과 결과를 tibble로 저장

한계효과는 플롯을 이용해서 보여줄 때, 가장 직관적으로 이해할 수 있습니다. 앞에서 구한 두 결과, 승산비로 계수값을 환산한 tidy_model_or과 한계효과를 구한 marginal_effects_df를 이용해서 플롯을 그려보겠습니다.

test_or <- tidy_model_or %>%
  dplyr::filter(term != "(Intercept)") %>% 
  select(variable = term, estimate = odds_ratio, ll, ul) %>% 
  mutate(Model = 'Odds Ratios')

test_ame <- marginal_effects_df %>%
  select(variable = factor, estimate = AME, ll = lower, ul = upper) %>%
  mutate(Model = 'Marginal Effects')

두 변수를 하나의 tibble로 결합하겠습니다.

data_for_graph <- bind_rows(test_or, test_ame)

그리고 이제 예측변수의 증가에 대한 종속변수 변화를 보여주는 승산비 모델과 한계효과 모델을 하나의 플롯으로 보여주되, 서로 그 기준을 바꾸어서 보여줄 것입니다.

  • 첫 번째 플롯은 승산비의 단위에 맞추어서 \(y\)축을 설정한 것이고, 두 번째 플롯은 한계효과의 단위에 맞춘 것입니다.

  • 여기서는 예제라서 두 모델을 하나의 플롯으로 보여주었지만, 만약 논문에 담을 것이라면 Logit 분석결과는 표로, 한계효과는 그래프로 보여주어 최대한 정보전달이 가능한 방식으로 구성할 수 있을 것입니다.

library(futurevisions)
data_for_graph %>%
  ggplot(aes(reorder(variable, estimate), estimate, color = Model)) +
  geom_point(position = position_dodge(width = .9)) +
  geom_linerange(aes(ymin = ll, ymax = ul), position = position_dodge(.9)) +
  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed') +
  labs(x = "Predictors", y = 'Estimate from Logit Model') +
  theme_bw() +
  scale_color_manual(values=futurevisions("mars")) 

## 잘 안보이니 한계효과를 따로 그려보겠습니다.

data_for_graph %>%
  ggplot(aes(reorder(variable, estimate), estimate, color = Model)) +
  geom_point(position = position_dodge(width = .9)) +
  geom_linerange(aes(ymin = ll, ymax = ul), position = position_dodge(.9)) +
  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed') +
  labs(x = "Predictors", y = 'Estimate from Logit Model') +
  theme_bw() + 
  scale_color_manual(values=futurevisions("mars")) + 
  coord_cartesian(ylim = c(0, .25))

R에서 Probit 모델 추정하기

Probit 모델은 Logit 모델과 유사하게 종속변수가 이항변수일 경우에 사용할 수 있으며, 당연히 그 결과도 Logit 모델과 매우 유사합니다.

  • 다만 종속변수가 이항변수일 경우, 대개 Logit 모델이 Probit 모델보다 선호되는데, 그 이유는 Logit 모델의 통계치들이 이용하기 편하기 때문입니다.

  • 또한 Logit 분석의 계수값은 지수화할 경우 승산비로 나타내어 해석할 수 있는 반면에, Probit 계수는 그렇게 해석할 수가 없습니다.

    • Logit 은 자연로그 분포를 이용하여 추정하는 방법인 반면에Probit은 정규분포의 CDF의 적분을 이용하는 추정방법이기 때문입니다.

  • Logit 모델과 Probit 모델이 서로 각자에게 우위를 갖는 (매우 미세하게) 특정한 사례들이 있을 수는 있지만, 대개의 경우 두 모델은 매우 유사한 결과를 보여주므로 대세를 따라 Logit으로 추정하는 것을 권합니다.

그럼 아까 Logit 모델 추정을 위해 사용했던 동일한 데이터를 Probit 모델 용으로 불러와보겠습니다.

df.probit <- pick("heart_data.dta")

## Probit 모델은 binomial(link = '')에서 logit이 probit으로 바뀐다는 것뿐입니다.
probit_model <- glm(
  death ~ anterior + hcabg + kk2 + kk3 + kk4 + 
    age2 + age3 + age4, 
  data = df.probit,
  family = binomial(link = 'probit')
)

summary(model) # Logit model
## 
## Call:
## glm(formula = death ~ anterior + hcabg + kk2 + kk3 + kk4 + age2 + 
##     age3 + age4, family = binomial(link = "logit"), data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3670  -0.3276  -0.1986  -0.1443   3.1807  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.0521     0.2586 -19.536  < 2e-16 ***
## anterior      0.6426     0.1676   3.835 0.000126 ***
## hcabg         0.7444     0.3530   2.109 0.034929 *  
## kk2           0.8117     0.1805   4.497 6.90e-06 ***
## kk3           0.7757     0.2691   2.883 0.003939 ** 
## kk4           2.6597     0.3560   7.471 7.96e-14 ***
## age2          0.4930     0.3102   1.589 0.111961    
## age3          1.5112     0.2662   5.676 1.38e-08 ***
## age4          2.1853     0.2718   8.039 9.06e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1484.6  on 4482  degrees of freedom
## Residual deviance: 1273.3  on 4474  degrees of freedom
##   (905 observations deleted due to missingness)
## AIC: 1291.3
## 
## Number of Fisher Scoring iterations: 7
summary(probit_model) # Probit model
## 
## Call:
## glm(formula = death ~ anterior + hcabg + kk2 + kk3 + kk4 + age2 + 
##     age3 + age4, family = binomial(link = "probit"), data = df.probit)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2888  -0.3285  -0.1932  -0.1307   3.2498  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.56975    0.10479 -24.522  < 2e-16 ***
## anterior     0.30029    0.07681   3.909 9.25e-05 ***
## hcabg        0.34890    0.17686   1.973  0.04852 *  
## kk2          0.37573    0.08536   4.402 1.07e-05 ***
## kk3          0.41049    0.13178   3.115  0.00184 ** 
## kk4          1.42881    0.20113   7.104 1.21e-12 ***
## age2         0.18332    0.12654   1.449  0.14741    
## age3         0.64877    0.11033   5.880 4.10e-09 ***
## age4         1.00223    0.11685   8.577  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1484.6  on 4482  degrees of freedom
## Residual deviance: 1268.6  on 4474  degrees of freedom
##   (905 observations deleted due to missingness)
## AIC: 1286.6
## 
## Number of Fisher Scoring iterations: 7
## 두 계수값이 유사한지 잘 확인하기 어렵다면 일단 데이터를 tidy하게 정리합니다.
tidy_probit <- probit_model %>% tidy(); tidy_probit
## 마찬가지로 신뢰구간도 구하고
tidy_probit <- tidy_probit %>%
  bind_cols(
    probit_model %>% broom::confint_tidy() %>%
      select(ll = conf.low, ul = conf.high)
  )
tidy_probit
## 변수 이름에 '.'을 '_'로 모두 변경. 
tidy_probit <- tidy_probit %>%
  rename_at(vars(contains('.')), function(x) str_replace(x, '\\.', '_'))

tidy_model # logit
tidy_probit # probit

Probit 모델도 예측확률을 구할 수 있습니다.

probit_new_data <- df %>% select(anterior, hcabg, kk2:kk4, age2:age4)
probit_predictions_r <- predict(
  probit_model, newdata = probit_new_data, type = 'response' 
  ## type = 'response' 옵션은 예측확률을 구하라는 옵션
  ## 이거 설정 안하면 그냥 한계효과만 죽 구합니다.
) %>% unname()

## 앞에서 미리 계산해두었던 Logit 모델의 예측확률과 Probit 모델의 예측확률을
## 하나의 tibble로 데이터로 만들어 줍니다.
model_predictions <- tibble(
  logit = predictions_r,
  probit = probit_predictions_r,
) %>% na.omit()

## tidyr 패키지의gather() 함수를 이용해서 Model이라는 변수로 long-shape 변환,
## 그리고 pred라는 변수에 값(value)들을 마찬가지로 long-shape로 변환합니다.
## Model 변수는 요인형으로 구분해줍니다.

model_predictions <- model_predictions %>%
  tidyr::gather(key = 'Model', value = 'pred') %>%
  mutate(
    Model = parse_factor(Model, levels = c('logit', 'probit'), include_na = F)
  )

이렇게 만든 model_predictions는 Logit과 Probit 모델 각각을 통해 구한 예측확률을 나타냅니다. group_by() 함수를 이용해 모델별로 pred를 정렬한 뒤 관측치 개수대로 목록(index) 변수를 만들고 ungroup으로 group_by()인한 tbl_group 상태를 해제하여 줍니다.

model_predictions <- model_predictions %>%
  group_by(Model) %>% 
  arrange(pred) %>% 
  mutate(index = 1:n()) %>% 
  ungroup()

그리고 \(x\)축을 관측치 넘버링, \(y\)축을 모델의 예측확률, 그리고 두 모델을 서로 다른 색으로 구분하여 주는 옵션을 통해 플롯을 그려주면 거의 유사한 패턴을 보임을 알 수 있습니다.

model_predictions %>%
  ggplot(aes(x = index, y = pred, color = Model)) +
    geom_path() +
    labs(y = 'Predictions from Model') +
    theme_bw() + 
  scale_color_manual(values=futurevisions("mars"))

이번에는 Probit model에서 한계효과를 구해 Logit 모델의 한계효과와 비교해보도록 하겠습니다. 먼저 Probit 모델의 평균한계효과와 그 신뢰구간을 구해서 tibble로 저장하고, Probit 모델의 한계효과임을 구분할 수 있게 Model 변수에 “Probit”이라고 입력하여 주었습니다.

probit_ame <- margins(probit_model)
summary(probit_ame)
probit_ame <- summary(probit_ame) %>%
  as.data.frame() %>% as_tibble()

probit_ame <- probit_ame %>%
  select(variable = factor, estimate = AME, ll = lower, ul = upper) %>%
  mutate(Model = 'Probit')

그리고 Logit 모델의 평균한계효과와 그 신뢰구간도 동일하게 구한 뒤, Model 변수를 만들어 “Logit”이라고 입력하여 주었습니다.

logit_ame <- test_ame # from logit model before
logit_ame <- logit_ame %>% mutate(Model = 'Logit')

logit_ame
probit_ame
## 당연히 두 tibble을 하나로 결합해 준다.
data_for_comparison_graph <- bind_rows(logit_ame, probit_ame)

두 모델의 평균한계효과를 비교해보았을 때, 추세와 통계적 유의성에 있어서 큰 차이가 나타나지 않는 것을 확인할 수 있습니다.

data_for_comparison_graph %>%
  ggplot(aes(reorder(variable, estimate), estimate, color = Model)) +
  geom_point(position = position_dodge(width = .9)) +
  geom_linerange(aes(ymin = ll, ymax = ul), position = position_dodge(.9)) +
  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed') +
  labs(x = "Predictors", y = 'Maringal Effect from Model') +
  theme_bw() +
  scale_color_manual(values=futurevisions("mars"))