1.1 Linear Model의 정의

• 반응변수와 독립변수간 관계를 나타낼 때, 독립변수와 모수간 선형결합의 형태를 띄는 것을 의미합니다.
• 또한, 반응변수가 정규분포 형태를 갖는 모델입니다.
• 예를 들어, Linear regression에서는 아래와 같은 식을 가지게됩니다.
  • \(Y = X^{T}\beta + \epsilon\)
  • \(Y\) : 반응변수. 연속적인 값
  • \(X\) : 1과 독립변수. matrix (아래 X 행렬의 모습을 확인하세요)
  • \(\beta\) : unknown parameters
  • 이 때, \(X^{T}\beta = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + ... + \beta_{p}X_{p}\)
  • \(\epsilon\) : 평균이 0인 random variable
• 여기서 Linear Model의 경우 Linear regression을 주로 이야기합니다.
• Time Series(ARIMA 모형 등)는 다루지 않습니다.

1.2 Linear regression model과 Least Squares

• 데이터가 아래와 같다고 가정
  • \((x_{1}, y_{1}), ... (x_{N}, y_{N})\) : N개의 training data
  • \(x_{i} = (x_{i1}, x_{i2}, ..., x_{ip})^{T}\) : i-th 관측치, p는 독립변수의 갯수
• 모델은 다음과 같음
  • \(f(X) = \beta_{0} + \Sigma_{j=1}^{p}X_{j}\beta_{j}\)
• 여러 방법이 존재하지만, 가장 대중적인 방법은 least squares estimation(최소제곱추정법) 입니다.
  • 아래 식을 최소화하는 \(\beta = (\beta_{0}, \beta_{1}, ..., \beta_{p})^{T}\)를 찾는 것입니다.

\[ RSS(\beta) = \Sigma_{i=1}^{N}(y_{i} - f(x_{i}))^2 \\ \quad = \Sigma_{i=1}^{N}(y_{i} - \beta_{0} - \Sigma_{j=1}^{p}x_{ij}\beta_{j})^2 \]

• 그렇다면, 어떻게 위 식을 최소화 할 수 있을까요?
  • \(X : N \times (p + 1)\)의 행렬을 정의, 이 때 1-벡터의 위치는 가장 첫번째에 위치합니다.
  • \(\mathbf y\) : 길이가 N인 1차원 벡터입니다.
  • \(\beta\) : 추정하고자하는 coefficient 벡터로 길이가 p + 1인 1차원벡터입니다.

\[ X = \begin{bmatrix} 1 & x_{11} & ... & x_{1p} \\ 1 & x_{21} & ... & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{N1} & ... & x_{Np} \end{bmatrix} \]

\[ \mathbf y = \begin{bmatrix} y_{1} \quad y_{2} \quad\cdots\quad y_{N} \end{bmatrix}^{T} \]

\[ \mathbf \beta = \begin{bmatrix} \beta_{0} \quad \beta_{1} \quad \beta_{2} \quad\cdots\quad \beta_{p} \end{bmatrix}^{T} \]

• \(RSS(\beta)\)를 행렬식으로 다시 표현

\[ RSS(\beta) = (\mathbf y - X\mathbf\beta)^{T}(\mathbf y - X\mathbf\beta) \]

• \(RSS(\beta)\)는 (p + 1)개의 parameters에 대한 이차형식(quadratic)입니다.

• 이를 \(\beta\) 관점에서 미분하면 다음과 같은 식을 얻을 수 있습니다.

\[ \frac{\partial RSS}{\partial \mathbf\beta} = -2X^{T}(\mathbf y - X\mathbf\beta)\\ \frac{\partial^{2} RSS}{\partial \mathbf\beta \partial \beta^{T}} = 2X^{T}X \]

• \(X\)가 Full column rank를 가지고, \(X^{T}X\)가 positive definite이면, 1차 미분을 0으로 셋팅할 수 있습니다.

\[ X^{T}(\mathbf y - X\mathbf\beta) = 0 \] - 우리의 목적은 \(\beta\)를 추정하는 것이므로 위 방정식에서 다음과 같이 unique solution을 얻을 수 있습니다.

\[ \hat{\mathbf \beta} = (X^{T}X)^{-1}X^{T}\mathbf y \space ... (a) \]

• \(\hat\beta\)가 training data를 이용해서 추정한 parameters가 되며, 이를 이용해 예측을 하려면 아래와 같이 계산을 해야합니다.

• 만일 test data를 \(X^{'}\)이라 하고 예측값을 \(\mathbf{\hat{y}}_{test}\)라 한다면,

\[ \begin{align} \mathbf{\hat{y}}_{test} = X^{'}\hat{\mathbf\beta} \end{align} \]

기하학적 의미

• 아래 그림은 설명변수가 2개일 때, least squares regression의 \(N\)차원의 기하학적 구조입니다.

• 노란색 평면은 \(x_{1}, x_{2}\)와 \(\beta_{0}, \beta_{1}, \beta_{2}\)의 선형결합으로 이루어져 있습니다.
  • 즉, 평면상에 있는 값이 바로 \(\hat{\mathbf y} ( = \hat\beta_{0} + \hat\beta_{1}x_{1} + \hat\beta_{2}x_{2} )\)입니다.
• 그렇다면 \(||\mathbf y - \hat{\mathbf y}||^2\)를 최소화하려면 \(\hat{\mathbf y}\)의 위치는 어떻게 되어야할까요?
  • \(\mathbf y - \hat{\mathbf y}\) 벡터와 노란색 평면이 수직이 되게끔 \(\beta_{0}, \beta_{1}, \beta_{2}\)를 구하면 됩니다.

한계점

• 위에서 \(X^{T}X\)의 역행렬을 이용해서 parameters를 추정하는데 만일 역행렬이 존재하지 않는 경우는 무엇이 있을까요?
  • 1. \(\mathbf X\)의 설명변수들이 Linearly Independent가 아닌 경우 (e.g. \(x_{2} = 3x_{1}\))
    • 즉, 두 변수간의 Correlation이 높을 때
  • 2. \(p > N\)인 경우

• 1)은 분석을 통해 Correlation이 높은 변수들을 제거하여 사용할 수도 있고, R의 lm 함수에서는 그러한 변수에 대해서는 자동적으로 제외합니다.

• 2)는 변수 선택 혹은 regularization을 수행하여 모델링을 합니다.

1.3 matrix 연산을 통한 모델링

# Toy Example을 이용한 linear regression 직접 구현
X <- matrix(rnorm(10000 * 5), ncol = 5)
y <- apply(X * matrix(c(1.3, 5.1, 3.4, 1.7, 4.2), nrow = nrow(X), ncol = 5, 
    byrow = TRUE), 1, sum) + 0.3

toy_df <- cbind(data.frame(X), y)
str(toy_df)

'data.frame':   10000 obs. of  6 variables:
 $ X1: num  0.0515 0.2118 0.6078 -0.604 -0.49 ...
 $ X2: num  -0.4881 0.5306 -2.5159 0.0533 -0.2573 ...
 $ X3: num  -0.448 -0.559 -0.71 0.984 -0.547 ...
 $ X4: num  1.131 0.355 -2.117 0.501 2.021 ...
 $ X5: num  0.0287 -1.5253 -2.5212 0.221 0.0264 ...
 $ y : num  -1.601 -4.4212 -28.3419 4.911 0.0371 ...

• Toy Example을 만들어서 간단하게 구현해봤습니다.

• 실제 정답은 \(\beta_{0} ~ \beta_{5}\)까지 순서대로 0.3, 1.3, 5.1, 3.4, 1.7, 4.2 입니다.

X_mat <- cbind(rep(1, nrow(X)), X)
beta_est <- solve(t(X_mat) %*% X_mat) %*% t(X_mat) %*% y
row.names(beta_est) <- paste0("b", 0:5)
beta_est

   [,1]
b0  0.3
b1  1.3
b2  5.1
b3  3.4
b4  1.7
b5  4.2

• Least Squares를 통해 추정된 \(\beta\)값을 구합니다.

• R의 lm함수를 이용하면 쉽게 Linear regression 모델을 생성할 수 있습니다.

1.4 R의 lm함수를 이용한 간단한 모델링

# lm함수를 이용한 linear regression 수행
toy_lm <- lm(y ~ ., data = toy_df)
coef(toy_lm)

(Intercept)          X1          X2          X3          X4          X5 
        0.3         1.3         5.1         3.4         1.7         4.2 

• 이번에는 R에서의 lm함수를 이용하여 \(\beta\)값을 추정해보았습니다.

summary(toy_lm)

Call:
lm(formula = y ~ ., data = toy_df)

Residuals:
       Min         1Q     Median         3Q        Max 
-3.969e-13 -2.100e-15 -8.000e-16  8.000e-16  4.736e-12 

Coefficients:
             Estimate Std. Error   t value Pr(>|t|)    
(Intercept) 3.000e-01  4.979e-16 6.026e+14   <2e-16 ***
X1          1.300e+00  4.999e-16 2.601e+15   <2e-16 ***
X2          5.100e+00  4.924e-16 1.036e+16   <2e-16 ***
X3          3.400e+00  4.968e-16 6.843e+15   <2e-16 ***
X4          1.700e+00  4.996e-16 3.403e+15   <2e-16 ***
X5          4.200e+00  4.960e-16 8.468e+15   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.978e-14 on 9994 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 4.882e+31 on 5 and 9994 DF,  p-value: < 2.2e-16

• summary함수를 통해 Linear regression model의 결과를 더 자세히 확인할 수 있습니다.

2.1 Generalized Linear Model의 정의

• 지금까지 Linear Model을 살펴봤습니다.

• 이번에는 Generalized Linear Model(이하 GLM)에 대해서 이야기를 해보겠습니다.

• GLM은 말 그대로 일반화된 선형 모델인데,

  • 일반화된 이라는 의미는 (제가 이해하기로는) 반응변수의 분포가 정규분포 이외에도 다른 분포들을 갖는 경우를 한번에 표현하기 위한 형용어구입니다. (이 부분은 피드백 부탁드리겠습니다.)

  • 이 때, 반응 변수의 분포가정과 Link function에 따라서 GLM의 성능이 달라질 수 있습니다.
    • 분포가정은 흔히 지수족 분포로 합니다. (Binomial, Poisson, multinomial, normal, …)
    • Link function의 경우 분포 가정에 맞추어서 설정합니다. (Logit, Identity, Log, …)

  • 예를 들어 반응 변수가 0 or 1, Good or Bad 등 Binary한 값들로 구성되어있는데 생뚱맞게 Gamma 분포로 가정한다면 결과가 이상할 수 있습니다.

  • Link function의 경우, 반응변수와 ’독립변수와 모수간 선형결합’간의 관계를 가정할 때 사용합니다.
    • 즉, \(E(Y) = \mu\) 이라 할 때,

\[g(\mu) = X\beta\]

• 이 때, \(g\) 함수가 Link function입니다.

• 그렇다면, Linear regression은 GLM관점에서 분포가정과 Link function이 각각 어떻게 셋팅되어야 할까요?
  • 반응변수는 정규분포라 가정하고,
  • Link function은 identity function을 가지면 됩니다.
• 우리가 주로 Binary Classification에서 사용하는 Logistic regression의 경우,
  • 반응변수는 Binomial 분포라 가정하고,
  • Link function은 Logit function 이라 하면 됩니다.

2.2 Generalized Linear model과 Maximum Likelihood estimation

• GLM에서의 parameters를 추정하는 방법은 여러가지이지만 보통은 Maximum Likelihood를 이용합니다.

• 그중에서도 주로 사용되는 Logistic regression을 MLE로 풀어보도록 하겠습니다.

  • \(X : N \times (p + 1)\)의 행렬을 정의, 이 때 1-벡터의 위치는 가장 첫번째에 위치합니다.
  • \(\mathbf y\) : 길이가 N인 1차원 벡터입니다. 개별 값들은 0 또는 1을 가집니다.
  • \(\beta\) : 추정하고자하는 coefficient 벡터로 길이가 p + 1인 1차원 벡터입니다.
  • \(g\) : Link function이며, 여기서는 Logit 함수입니다.

\[ X = \begin{bmatrix} 1 & x_{11} & ... & x_{1p} \\ 1 & x_{21} & ... & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{N1} & ... & x_{Np} \end{bmatrix} \]

\[ \mathbf y = \begin{bmatrix} y_{1} \quad y_{2} \quad\cdots\quad y_{N} \end{bmatrix}^{T} \]

\[ \mathbf \beta = \begin{bmatrix} \beta_{0} \quad \beta_{1} \quad \beta_{2} \quad\cdots\quad \beta_{p} \end{bmatrix}^{T} \]

\[ g(p) = ln( \frac{p}{1-p} ) \]

• MLE를 구하기 위해서는 Likelihood function이 필요한데, 이는 probability density function(이하 pdf)을 parameters관점에서 바라본 함수입니다.

• \(y_{i}, (i = 1, ..., N)\)는 반응변수의 i번째 관측치로써, \(Binomial(1, p)\)를 따른다고 가정합니다.

• 그러면 \(y_{i}\)의 pdf는 다음과 같습니다.

\[ P(y_{i}|p) = p^{y_{i}}(1-p)^{y_{i}} \] - 이 때 N개의 관측치들의 pdf는 다음과 같습니다.

\[ \begin{align} P(\mathbf y|p) & = P(y_{1}, y_{2}, ..., y_{N}|p) \\ & = P(y_{1}|p) P(y_{2}|p) ... P(y_{N}|p) \\ & = \Pi_{i = 1}^{N}P(y_{i}|p) \\ & = \Pi_{i = 1}^{N}p^{y_{i}}(1-p)^{1-y_{i}} \\ & = \Pi_{i = 1}^{N} (\frac{p}{1-p})^{1-y_{i}} (1-p) \end{align} \]

• regression은 \(y_{i}\)가 아니라 \(y_{i}\)의 평균값을 예측하는 것이 목적입니다.

• \(y_{i}\)의 분포를 \(Binomial(1, p)\)라고 가정하였으므로,
  • \(E(y_{i}) = p\)에 대해서 모델링을 수행하는 것입니다.

• 그리고 Logit 함수를 Link function으로 이용하므로, 아래와 같은 식을 이용하여 모델링을 합니다.

\[ \begin{align} g(p) & = ln( \frac{p}{1-p} ) \\ & = X\beta \end{align} \]

• 식을 조금만 더 변형시켜볼까요?

\[ \frac{p}{1-p} = exp(X\beta) \\ 1-p = \frac{1}{1 + exp(X\beta)} \]

• 위 식을 pdf에 대입을 하면 아래와 같습니다.
  • 드디어 Likelihood function을 얻었습니다!

\[ L(\beta|\mathbf y, X) = \Pi_{i=1}^{N} exp(X_{i}\beta)^{y_{i}} \frac{1}{1+ exp(X_{i}\beta)} \]

• \(X_{i}\)는 행렬 \(X\)의 i번째 행을 의미합니다.

• Likelihood 함수에 log를 취해보겠습니다.

\[ \begin{align} l(\beta|\mathbf y, X) & = log(L(\beta|\mathbf y, X)) \\ & = log(\Pi_{i=1}^{N} exp(X_{i}\beta)^{y_{i}} \frac{1}{1+ exp(X_{i}\beta)}) \\ & = \Sigma_{i = 1}^{N}y_{i}X_{i}\beta - \Sigma_{i = 1}^{N}log(1+exp(X_{i}\beta)) \end{align} \]

• 위와 같은 식을 Log-Likelihood function이라고 합니다.

• 위 식을 최대화하는 \(\beta\)를 찾는 것이 바로 Logistic regression에서 모델링을 수행하는 것입니다.

• 1차 미분과 2차 미분을 수행하여 찾을 수 있는데 이 때 사용되는 알고리즘들이 있습니다.
  • 1차 미분과 2차 미분은 Log-Likelihood 함수가 최대화되는 지점을 찾기 위해서 계산합니다.
  • 그런데 2차 미분은 계산량이 너무 많아서 다른 알고리즘들을 사용하여 근사적으로 구합니다.
    • 대표적인 방법으로 Iteratively Reweighted Least Squares Method(IRLSM)이 있습니다.
  • 제가 아직 IRLSM을 잘몰라서 본 문서에는 따로 기재하지 않았습니다.

2.3 R의 glm함수를 이용한 간단한 모델링

• 지금까지 GLM의 이론적인 부분을 간략하게(?) 살펴봤습니다.
• 너무 빙 돌아서 온 것 같지만, 그냥 사용하는 것보단 알고 사용하는 것이 삽질도 덜 할 수 있습니다.
• 이번에는 R의 glm함수를 이용해서 Logistic regression 모델을 만들어보겠습니다.

Load example data

# 예제 데이터 가져오기
library(ISLR)
library(dplyr)
data("Default", package = "ISLR")
glimpse(Default)

Observations: 10,000
Variables: 4
$ default <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, No...
$ student <fct> No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, N...
$ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 919...
$ income  <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496,...

• ISLR 패키지에서 Defalut 데이터를 가져왔습니다.

• Default 데이터는 고객들의 신용카드 채무이행/불이행 관련 데이터로 아래와 같은 변수들을 가지고 있습니다.
  • 10,000 X 4의 차원을 가지는 데이터
  • default (factor) : 반응변수. 채무불이행관련 변수(Yes : 채무불이행 / No : 채무이행)
  • student (factor) : 독립변수. 학생인지 여부
  • balance (numeric) : 독립변수. 월별 결제 후 고객이 신용카드에 남아있는 평균 잔액
  • income (numeric) : 독립변수. 고객의 소득

GLM을 이용한 모델링

glm_default <- glm(default ~ ., data = Default, family = "gaussian")

• 위 코드를 실행하면 아래와 같은 에러가 뜰겁니다.

Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, : 'y' 내에 NA/NaN/Inf가 있습니다

• default는 factor형 자료구조를 가지는데 gaussian으로 적용하려하니 일어나는 일입니다.

• family = binomial로 설정해서 수행하면 Logistic regression 모델을 만들 수 있습니다.

glm_default <- glm(default ~ ., data = Default, family = "binomial")
summary(glm_default)

Call:
glm(formula = default ~ ., family = "binomial", data = Default)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4691  -0.1418  -0.0557  -0.0203   3.7383  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.087e+01  4.923e-01 -22.080  < 2e-16 ***
studentYes  -6.468e-01  2.363e-01  -2.738  0.00619 ** 
balance      5.737e-03  2.319e-04  24.738  < 2e-16 ***
income       3.033e-06  8.203e-06   0.370  0.71152    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2920.6  on 9999  degrees of freedom
Residual deviance: 1571.5  on 9996  degrees of freedom
AIC: 1579.5

Number of Fisher Scoring iterations: 8

• Logistic regression 모델을 만들고 summary를 통해서 모델의 정보도 확인이 가능합니다.

• 그렇다면, 이러한 것을 h2o에서도 할 수 있을까요?

2.4 h2o를 이용한 간단한 모델링

• 이번에는 h2o를 이용하여 Logistic regression 모델을 만들어보도록 하겠습니다.

library(h2o)
h2o.init(nthreads = 10, max_mem_size = "16g")

H2O is not running yet, starting it now...

Note:  In case of errors look at the following log files:
    C:\Users\AGILES~1\AppData\Local\Temp\RtmpqgP6E8/h2o_agilesoda_started_from_r.out
    C:\Users\AGILES~1\AppData\Local\Temp\RtmpqgP6E8/h2o_agilesoda_started_from_r.err


Starting H2O JVM and connecting: . Connection successful!

R is connected to the H2O cluster: 
    H2O cluster uptime:         8 seconds 313 milliseconds 
    H2O cluster timezone:       Asia/Seoul 
    H2O data parsing timezone:  UTC 
    H2O cluster version:        3.20.0.2 
    H2O cluster version age:    13 days  
    H2O cluster name:           H2O_started_from_R_agilesoda_myv006 
    H2O cluster total nodes:    1 
    H2O cluster total memory:   14.19 GB 
    H2O cluster total cores:    12 
    H2O cluster allowed cores:  10 
    H2O cluster healthy:        TRUE 
    H2O Connection ip:          localhost 
    H2O Connection port:        54321 
    H2O Connection proxy:       NA 
    H2O Internal Security:      FALSE 
    H2O API Extensions:         Algos, AutoML, Core V3, Core V4 
    R Version:                  R version 3.4.4 (2018-03-15) 

y <- "default"
x <- names(Default)[!(names(Default) %in% y)]

Default.h2o <- as.h2o(Default, destination_frame = "Default.h2o")

  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

h2o_data <- h2o.splitFrame(Default.h2o, ratios = 0.7, destination_frames = c("train_dt.h2o", 
    "test_dt.h2o"))

• List 형태로 나오며, 첫번쨰가 training set, 두번째가 test set입니다.

glm_default_h2o <- h2o.glm(x = x, y = y, training_frame = h2o_data[[1]], family = "binomial", 
    model_id = "glm_default_h2o")

  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |=================================================================| 100%

glm_default_h2o

Model Details:
==============

H2OBinomialModel: glm
Model ID:  glm_default_h2o 
GLM Model: summary
    family  link                                regularization
1 binomial logit Elastic Net (alpha = 0.5, lambda = 1.265E-4 )
  number_of_predictors_total number_of_active_predictors
1                          4                           4
  number_of_iterations training_frame
1                    6   train_dt.h2o

Coefficients: glm coefficients
        names coefficients standardized_coefficients
1   Intercept   -11.435978                 -6.357521
2  student.No     0.320279                  0.320279
3 student.Yes    -0.301216                 -0.301216
4     balance     0.005840                  2.826461
5      income     0.000007                  0.088041

H2OBinomialMetrics: glm
** Reported on training data. **

MSE:  0.0213689
RMSE:  0.146181
LogLoss:  0.07756588
Mean Per-Class Error:  0.2633005
AUC:  0.9533673
Gini:  0.9067346
R^2:  0.338751
Residual Deviance:  1081.113
AIC:  1091.113

Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
         No Yes    Error       Rate
No     6658  78 0.011580   =78/6736
Yes     120 113 0.515021   =120/233
Totals 6778 191 0.028412  =198/6969

Maximum Metrics: Maximum metrics at their respective thresholds
                        metric threshold    value idx
1                       max f1  0.334653 0.533019 117
2                       max f2  0.148635 0.597946 203
3                 max f0point5  0.448191 0.591241  82
4                 max accuracy  0.448191 0.973597  82
5                max precision  0.980341 1.000000   0
6                   max recall  0.000939 1.000000 390
7              max specificity  0.980341 1.000000   0
8             max absolute_mcc  0.334653 0.521241 117
9   max min_per_class_accuracy  0.040280 0.881384 296
10 max mean_per_class_accuracy  0.037306 0.888292 300

Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`