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

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

library(dmtr)
library(dplyr)

데이터 {#multiple-linear-regression-data}

data(biometric, package = "dmtr")

knitr::kable(
  biometric, 
  booktabs = TRUE,
  align = c('r', 'r', 'r'),
  caption = '나이, 키, 몸무게 데이터'
)

다중회귀모형 추정 {#multiple-linear-regression}

아래와 같이 $n$개의 객체와 $k$개의 독립변수($\mathbf{x}$)로 이루어지고 하나의 종속변수($y$)로 이루어진 선형 회귀모형을 정의하자.

\begin{equation} y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik} + \epsilon_i (#eq:multiple-linear-regression) \end{equation}

이 때, 오차항 $\epsilon_i$은 서로 독립이고 동일한 정규분포 $N(0, \sigma^2)$을 따른다.

위 회귀모형은 아래와 같이 행렬의 연산으로 표한할 수 있다.

\begin{equation} \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} (#eq:multiple-linear-regression-matrix) \end{equation}

이 때,

[ \mathbf{y} = \left[ \begin{array}{c} y_1 \ y_2 \ y_3 \ \vdots \ y_n \end{array} \right] ]

[ \mathbf{X} = \left[ \begin{array}{c c c c c} 1 & x_{11} & x_{12} & \cdots & x_{1k}\ 1 & x_{21} & x_{22} & \cdots & x_{2k}\ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & x_{n1} & x_{n2} & \cdots & x_{nk} \end{array} \right] ]

[ \boldsymbol{\beta} = \left[ \begin{array}{c} \beta_0 \ \beta_1 \ \beta_2 \ \vdots \ \beta_k \end{array} \right] ]

[ \boldsymbol{\epsilon} = \left[ \begin{array}{c} \epsilon_1 \ \epsilon_2 \ \epsilon_3 \ \vdots \ \epsilon_n \end{array} \right] ]

로 정의되며,

[ E[\boldsymbol{\epsilon}] = \mathbf{0}, \, Var[\boldsymbol{\epsilon}] = \sigma^2 \mathbf{I} ]

이다.

회귀계수추정 {#regression-coefficient-estimate}

다중회귀모형에서 회귀계수들의 추정은 최소자승법(least squares method)에 근거하고 있다. 최소화시킬 오차항의 제곱합 $Q$는 다음과 같이 표현된다.

\begin{equation} Q = \sum_{i = 1}^{n} \left(y_i - \left(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik}\right) \right) ^ 2 \end{equation}

최소자승법에 의한 회귀계수의 추정은 제곱합 $Q$를 각 $\beta_j$에 대하여 편미분하고 이를 0으로 하는 다음과 같은 연립방정식을 풀어 $\hat{\beta}_j$들을 구하는 것이다.

\begin{eqnarray} \frac{\partial Q}{\partial \beta_0} &=& - 2 \sum_{i = 1}^{n} 1 \left(y_i - \left(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik}\right) \right) = 0\ \frac{\partial Q}{\partial \beta_j} &=& - 2 \sum_{i = 1}^{n} x_{ij} \left(y_i - \left(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik}\right) \right) = 0, \, j = 1, \cdots, k \end{eqnarray}

여기에서

[ x_{i0} = 1, \, i = 1, \cdots, n ]

이라 하면, 위 오차항의 제곱합과 회귀계수에 대한 편미분식은 아래와 같이 정리할 수 있다.

\begin{equation} Q = \sum_{i = 1}^{n} \left(y_i - \left(\beta_0 x_{i0} + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik}\right) \right) ^ 2 (#eq:multiple-linear-regression-sse) \end{equation}

\begin{equation} \frac{\partial Q}{\partial \beta_j} = - 2 \sum_{i = 1}^{n} x_{ij} \left(y_i - \left(\beta_0 x_{i0} + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik}\right) \right), \, j = 0, \cdots, k (#eq:multiple-linear-regression-gradient) \end{equation}

식 \@ref(eq:multiple-linear-regression-gradient)는 아래와 같이 행렬식으로 표현할 수 있다.

[ \mathbf{X}^\top \left( \mathbf{y} - \mathbf{X}\mathbf{\beta} \right) ]

이를 다시 정리하면 아래와 같이 회귀계수를 추정할 수 있다.

\begin{eqnarray} - 2 \mathbf{X}^\top \left( \mathbf{y} - \mathbf{X}\hat{\mathbf{\beta}} \right) &=& 0\ \mathbf{X}^\top \mathbf{y} &=& \mathbf{X}^\top \mathbf{X} \hat{\mathbf{\beta}}\ \left(\mathbf{X}^\top \mathbf{X}\right) ^ {-1} \mathbf{X}^\top \mathbf{y} &=& \left(\mathbf{X}^\top \mathbf{X}\right) ^ {-1} \left(\mathbf{X}^\top \mathbf{X}\right) \hat{\mathbf{\beta}}\ &=& \hat{\mathbf{\beta}} \end{eqnarray}

행렬식을 이용하여, 주어진 예제 데이터에서 나이와 키로써 몸부게를 설명하는 회귀모형을 추정해보자.

y <- biometric %>% 
  pull(weight)

X <- biometric %>% 
  select(age, height) %>%
  mutate(`(Intercept)` = 1, .before = 1L) %>%
  as.matrix()

beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y

print(beta_hat)

추정된 회귀계수에서 각 편미분값이 0을 확인해보자.

t(X) %*% (y - X %*% beta_hat)

실제 얻어진 값은 정확히 0은 아니지만 0에 매우 근접하며, 정확히 0이 얻어지지 않는 이유는 컴퓨팅 측면의 문제로 볼 수 있다 (알고리즘, 소수점 표현 구조 등). 이러한 경우 dplyr::near() 함수를 이용하여 수치가 목표 수치에 충분히 근접하여 있는지를 확인할 수 있으며, 위 회귀계수 편미분값의 경우 모두 0에 충분히 근접함을 확인할 수 있다.

near(t(X) %*% (y - X %*% beta_hat), 0)

위와 같은 회귀계수 추정과정을 함수 fit_linear_regression()으로 구현하였으며, 결과 리스트값에서 betas 원소가 추정된 회귀계수를 나타낸다.

fit <- fit_linear_regression(biometric, weight, c(age, height))
fit$betas

위 추정된 회귀계수를 식 \@ref(eq:multiple-linear-regression-sse)에 대입하여 잔차제곱합(residual sum of squareds; $SSE$)을 계산해보자.

sse <- sum((y - (X %*% fit$beta)) ^ 2)
sse

이 때, 식 \@ref(eq:multiple-linear-regression)의 오차항 $\epsilon_i$의 분산 $\sigma ^ 2$의 추정은 위 잔차제곱합을 자유도 $(n - k - 1)$로 나눈 잔차평균제곱합(mean sequared error; $MSE$)을 이용한다.

[ \hat{\sigma} ^ 2 = MSE = \frac{SSE}{n - k - 1} ]

mse <- sse / (nrow(X) - ncol(X))
mse

함수 fit_linear_regression()의 결과 리스트값에서 mse 원소가 잔차평균제곱값을 나타내며, df 원소가 자유도 $(n - k - 1)$을 나타낸다.

fit$mse

fit$df

따라서, 이 두 결과값을 이용하여 잔차제곱합을 역으로 계산할 수 있다.

fit$mse * fit$df

모형에 따른 추론 {#multiple-linear-regression-inference}

회귀성 검정

다중회귀모형에 있어서 유의성 검정, 즉 회귀성 검정은 독립변수의 기울기에 해당하는 모든 회귀계수($\beta_0$은 제외)가 0의 값을 갖는가 또는 그렇지 않은가를 검정하는 것이다. 귀무가설과 대립가설은 아래와 같다.

\begin{eqnarray} H_0 &:& \beta_1 = \beta_2 = \cdots = \beta_k = 0\ H_1 &:& \beta_j \neq 0 \, \text{for at least one} \, j \in {1, 2, \cdots, k} \end{eqnarray}

독립변수의 기울기에 해당하는 회귀계수 중 적어도 하나가 0이 아니라고 판단되면 귀무가설이 기각되고 회귀성이 있다고 말하게 된다.

회귀성 검정을 위해, 앞 절에서 살펴본 잔차제곱합이며 $SSE$ 외에, 전체제곱합(total sum of squares; $SST$)와 회귀제곱합(regression sum of squares; $SSR$)을 아래와 같이 정의한다.

\begin{eqnarray} SST &=& \sum_{i = 1}^{n} \left(y_i - \bar{y}\right) ^ 2\ SSR &=& \sum_{i = 1}^{n} \left(\hat{y}_i - \bar{y}\right) ^ 2 \end{eqnarray}

이 때, $\bar{y} = \frac{1}{n} \sum_{i = 1}^{n} y_i$는 관측치 $y_i$의 평균을 나타낸다. $SST$, $SSR$, 그리고 $SSE$는 아래와 같은 관계를 지닌다.

\begin{eqnarray} SST &=& \sum_{i = 1}^{n} \left(y_i - \bar{y}\right) ^ 2\ &=& \sum_{i = 1}^{n} \left(y_i - \hat{y}{i} + \hat{y}{i} - \bar{y}\right) ^ 2\ &=& \sum_{i = 1}^{n} \left(\left(y_i - \hat{y}{i}\right) + \left(\hat{y}{i} - \bar{y}\right)\right) ^ 2\ &=& \sum_{i = 1}^{n} \left(y_i - \hat{y}{i}\right) ^ 2 + \sum{i = 1}^{n} \left(\hat{y}{i} - \bar{y}\right) ^ 2 + 2 \sum{i = 1}^{n} \left(y_i - \hat{y}{i}\right) \left(\hat{y}{i} - \bar{y}\right)\ &=& SSE + SSR + 2 \sum_{i = 1}^{n} \left(y_i - \hat{y}{i}\right) \left(\hat{y}{i} - \bar{y}\right)\ &=& SSE + SSR \end{eqnarray}

여기에서 회귀제곱합 $SSR$을 자유도 $k$로 나눈 값을 아래와 같이 회귀평균제곱($MSR$)으로 아래와 같이 정의하고,

\begin{equation} MSR = \frac{SSR}{k} \end{equation}

검정통계량으로 다음과 같은 $F$-값을 이용한다.

\begin{equation} F_0 = \frac{MSR}{MSE} \end{equation}

회귀성 검정의 귀무가설이 옳을 경우, 즉 모든 독립변수의 기울기에 해당하는 회귀계수가 0일 경우, 검정통계량 $F_0$는 다음과 같은 F-분포를 따른다.

\begin{equation} F_0 \sim F_{(k, n - k - 1)}, \, \text{if } \beta_1 = \beta_2 = \cdots = \beta_k = 0 \end{equation}

따라서, 다음과 같을 때 유의수준 $\alpha$에서 가설 $H_0$을 기각한다.

[ F_0 > F_{(\alpha, k, n - k - 1)} ]

위 예제 데이터에서 추정된 회귀모형에 대해 회귀성 검정을 수행해보자.

sst <- sum((y - mean(y)) ^ 2)
ssr <- sst - fit$mse * fit$df
msr <- ssr / (ncol(X) - 1)
f0 <- msr / fit$mse
alpha <- 0.05
f_alpha <- qf(alpha, (ncol(X) - 1), fit$df, lower.tail = FALSE)
f0 > f_alpha

위 결과에서 귀무가설은 기각되며, 유의수준 r alpha에서 유의한 모형이라 할 수 있다. 위 결과를 함수 anova_linear_regression()을 이용하여 분산분석표로 나타낼 수 있다.

anova_linear_regression(fit)

개별 회귀계수에 대한 $t$-검정

회귀계수벡터 $\boldsymbol{\beta}$의 추정량은 다음과 같은 기대치와 분산-공분산행렬을 갖는 다변량 정규분포를 따른다.

\begin{eqnarray} E\left[\hat{\boldsymbol{\beta}}\right] &=& \boldsymbol{\beta}\ Var\left[\hat{\boldsymbol{\beta}}\right] &=& \sigma^2 \left(\mathbf{X}^\top \mathbf{X}\right)^{-1} \end{eqnarray}

위 예제 데이터에 대하여 회귀계수벡터 추정량의 분산-공분산행렬을 구해보자.

fit$mse * solve(t(X) %*% X)

함수 fit_linear_regression()의 결과 리스트값에서 hessian 원소는 $2 \mathbf{X}^\top \mathbf{X}$ 값을 지니며, 따라서 hessian 원소를 2로 나누어 $\mathbf{X}^\top \mathbf{X}$에 대입할 때 위와 동일한 분산-공분산행렬을 보인다.

fit$mse * solve(fit$hessian / 2)

참고로, hessian 행렬은 식 \@ref(eq:multiple-linear-regression-sse)의 오차항 제곱합 $Q$를 최적화하는 과정에서 얻어지며, 오차항 제곱합 $Q$의 식을 변형하면 hessian 행렬 또한 달라진다. 예를 들어, $Q$ 대신 $\frac{1}{2}Q$를 최적화하게 되면, 최적해 $\hat{\boldsymbol{\beta}}$는 동일하나 hessian 행렬은 $2 \mathbf{X}^\top \mathbf{X}$가 아닌 $\mathbf{X}^\top \mathbf{X}$로 얻어진다.

위 분산-공분산행렬에서 대각원소들은 각 회귀계수의 분산을 나타내며, 따라서 각 대각원소의 제곱근이 각 회귀계수의 표준오차(standard error; $se\left(\hat{\beta}_j\right)$)를 나타낸다.

sqrt(diag(fit$mse * solve(fit$hessian / 2)))

이 표본오차는 fit_linear_regression()의 결과 리스트의 se 원소에 저장되어 있다.

fit$se

이 때, $\hat{\beta}_j$의 분포는 다음과 같다.

[ \frac{\hat{\beta}j - \beta_j}{se\left(\hat{\beta}_j\right)} \sim t{(n - k - 1)} ]

그러므로 $\beta_j$에 대한 $100(1 - \alpha)\%$ 신뢰구간은 아래와 같다.

[ \hat{\beta}j \pm se\left(\hat{\beta}_j\right)t{(\alpha/2, \, n - k - 1)} ]

각각의 회귀계수($\beta_0$ 포함)가 0의 값을 갖는가 또는 그렇지 않은가를 검정하기 위한 귀무가설과 대립가설은 아래와 같다.

\begin{eqnarray} H_0 &:& \beta_j = 0\ H_1 &:& \beta_j \neq 0 \end{eqnarray}

위 검정을 위해 검정통계량 $t$-값을 이용한다.

[ T_j = \frac{\hat{\beta}_j}{se\left(\hat{\beta}_j\right)} ]

위 통계량은 $\beta_j = 0$일 때 $t_{(n - k - 1)}$ 분포를 따르기 때문에, 유의수준 $0 < \alpha < 1$에서

[ \left| T_j \right| > t_{(2 / \alpha, \, n - k - 1)} ]

이면 귀무가설 $H_0$를 기각하게 된다.

위 예제 데이터에서 추정된 회귀모형의 회귀계수 각각에 대해 검정을 수행해보자.

alpha <- 0.05
t_stat <- fit$betas / fit$se
t_alpha <- qt(alpha / 2, fit$df, lower.tail = FALSE)
abs(t_stat) > t_alpha

위 결과에서 모든 각각의 회귀계수에 대해 귀무가설은 기각되며, 유의수준 r alpha에서 각 독립변수는 종속변수인 몸무게에 유의한 변수라는 것을 의미한다. 위 결과를 함수 ttest_linear_regression()을 이용하여 데이터 프레임으로 나타낼 수 있다.

ttest_linear_regression(fit)

개별 기울기 회귀계수에 대한 $F$-검정

$r$개의 독립변수 $x_1, x_2, \cdots, x_r$이 포함된 다중회귀모형에서 변수 $x_j$에 대한 Type $\text{II}$ 제곱합이란, $x_j$를 제외한 타 변수가 이미 포함된 모형에 $x_j$가 추가로 도입될 때 증가되는 회귀제곱합을 의미한다.

[ \Delta SSR(x_j \, | \, x_1, \cdots, x_{j - 1}, x_{j + 1}, \cdots, x_{r}) = SSR(x_1, \cdots, x_{r}) - SSR(x_1, \cdots, x_{j - 1}, x_{j + 1}, \cdots, x_{r}) ]

여기서 $SSR(x_1, \cdots, x_{r})$은 변수 $x_1, \cdots, x_{r}$이 포함된 모형의 회귀제곱합($SSR$)을 말한다.

위 예제에 대한 완전회귀모형에서 독립변수 age에 대한 Type $\text{II}$ 제곱합을 계산해보자.

ssr <- fit$sst * fit$rsq
fit_without_age <- fit_linear_regression(biometric, weight, height)
ssr_without_age <- fit_without_age$sst * fit_without_age$rsq
type2_age <- ssr - ssr_without_age
type2_age

Type $\text{II}$ 제곱합의 유의성 검정은 $F$ 검정을 통하여 이루어진다. 변수 $x_j$에 대한 $F$-값은 다음과 같이 산출된다.

[ \Delta F_j = \frac{\Delta SSR(x_j \, | \, x_1, \cdots, x_{j - 1}, x_{j + 1}, \cdots, x_{r})}{MSE(x_1, \cdots, x_{r})} ]

여기에서 $MSE(x_1, \cdots, x_{r})$은 변수 $x_1, \cdots, x_{r}$이 포함된 모형의 잔차평균제곱($MSE$)을 말한다.

변수 $x_1, \cdots, x_{r}$이 포함된 모형에 대해 아래와 같은 회귀계수 검정을 수행한다 하자.

\begin{eqnarray} H_0 &:& \beta_j = 0\ H_1 &:& \beta_j \neq 0 \end{eqnarray}

귀무가설 $H_0$가 참일 때, $\Delta F_j$는 아래와 같은 $F$ 분포를 따른다.

[ \Delta F_j \sim F_{(1, n - r + 1)} ]

따라서, 각 변수의 회귀계수에 대한 유의수준 $\alpha$에서의 $F$ 검정은, $\Delta F_j > F(\alpha, 1, n - r + 1)$일 때 귀무가설을 기각하여 변수 $x_j$가 유의한 설명력을 지닌다고 판단한다. 위 예에서 변수 age에 대한 $F$-값을 구한 뒤, $F$ 검정을 수행해보자.

alpha <- 0.05
f_age <- type2_age / fit$mse
f_age > qf(alpha, 1, fit$df, lower.tail = FALSE)

위 결과에서, 유의수준 r alpha에서 귀무가설은 기각되며, 따라서 height 변수가 이미 존재하는 모형에 age 변수를 추가할 때 종속변수를 추가로 유의하게 설명한다 할 수 있다.

이와 같은 $F$ 검정을 회귀모형의 모든 각각의 변수에 수행하는 함수 test_type2_linear_regression()을 아래와 같이 사용할 수 있다.

test_type2_linear_regression(biometric, weight, c(age, height))

변수 선택

종속변수를 설명하기 위한 독립변수들의 후보가 매우 많은 경우, 어떤 독립변수들의 조합에 의한 회귀모형이 가장 좋은가를 판단한다.

모형 성능 척도

고려하고 있는 회귀모형 또는 추정된 회귀식이 얼마나 데이터를 잘 반영하고 있는가를 알기 위해, 회귀모형의 결정계수(coefficient of determination) $R ^ 2$를 아래와 같이 정의한다.

[ R ^ 2 = \frac{SSR}{SST} = 1 - \frac{SSE}{SST} ]

즉, $R ^ 2$값은 전체제곱합 중 모형이 설명하는 제곱합의 비율로 해석되며, 0에서 1 사이의 값을 갖는다. 위 예제 데이터에서 추정된 회귀모형의 결정계수 값을 계산해보자.

sst <- sum((y - mean(y)) ^ 2)
1 - (fit$mse * fit$df) / sst

이 식에서, $SSE$는 독립변수의 개수가 증가함에 따라 감소하는 속성을 지니며, $SST$는 독립변수의 개수와 상관없이 일정하므로, $R ^ 2$는 독립변수 개수가 증가함에 따라 증가하는 속성을 지닌다. $R ^ 2$에 기반한 모형 평가는 종종 과적합(overfitting)의 문제를 수반하므로, 독립변수 개수에 따른 $R ^ 2$의 증가분에 대한 조정을 적용하는 수정결정계수(adjusted coefficient of determination) $R_{adj} ^ 2$가 아래와 같이 정의된다.

[ R_{adj} ^ 2 = 1 - \frac{SSE / (n - k - 1)}{SST / (n - 1)} ]

여기에서 분자 $SSE / (n - k - 1)$은 앞 절에서 살펴본 잔차평균제곱합 $MSE$ 이며, 분모는 $SST$를 그 자유도 $n - 1$으로 나눈 값이다.

1 - fit$mse / (sst / (length(y) - 1))

위 결정계수 $R ^ 2$ 및 수정결정계수 $R_{adj} ^ 2$는 아래와 같이 fit_linear_regression() 결과 리스트의 rsq 및 rsqadj 원소에 각각 저장된다.

print(fit$rsq)
print(fit$rsqadj)

부가적인 모형의 성능척도로, 맬로우즈(Mallows)가 개발한 $C_p$ 통계량이 종종 사용된다. 이는 $p$항 회귀모형이 완전모형에 얼마나 가까운지를 나타내는 척도로서, 다음과 같이 정의된다.

[ C_p = \frac{SSE_p}{MSE} - n + 2p ]

여기에서 $MSE$는 $k + 1$개의 회귀계수(($k$개의 기울기 계수와 1개의 절편)를 이용한 완전모형에서 얻어지는 $\sigma^2$의 추정량이며, $SSE_p$는 $(p - 1)$개의 독립변수에 대한 회귀계수와 1개의 절편으로 이루어진 모형을 추정할 때 얻어지는 잔차제곱합이다 ($p - 1 < k$).

예를 들어, 위 예제 데이터에서 나이를 제외하고 키 하나의 변수만을 이용하여 몸무게를 설명하는 회귀모형을 아래와 같이 추정해보자.

fit_reduced <- fit_linear_regression(biometric, weight, height)

이 때, $C_p$값은 아래와 같이 얻어진다.

(fit_reduced$mse * fit_reduced$df) / fit$mse - 
  fit_reduced$n + 2 * (fit_reduced$n - fit_reduced$df)

참고로, $C_p$는 다음과 같은 예측치 평균오차제곱의 추정량이다.

[ \Gamma_p = \frac{1}{\sigma^2} \sum_{i = i}^{n} E\left[\left(\hat{y}_i - E\left[y_i\right]\right)^2\right] ]

완전모형($p = k + 1$)에 대한 $C_p$값은 $k + 1 = p$이다. 또한, $p$항 회귀모형이 완전모형에 비해 편의(bias)가 없다면, $C_p$의 기대값은 $p$가 된다.

\begin{eqnarray} E\left[C_p \, | \, bias = 0\right] &\approx& \frac{E\left[SSE_p \, | \, bias = 0\right]}{\sigma^2} - n + 2p\ &=& \frac{(n - p)\sigma^2}{\sigma^2} - n + 2p\ &=& p \end{eqnarray}

$C_p$값이 $p$에 가까울수록 편의가 적은 모형이며, 편의가 적은 모형인 동시에 변수의 수가 적은 모형이 가장 적절한 모형이라 할 수 있다.

모든 가능한 조합의 회귀분석

이 방법은 모든 가능한 독립변수들의 조합에 대한 회귀모형을 분석하여 가장 적합한 회귀모형을 선택하는 방법이다.

variables <- c("age", "height")
variables_in_model <- vector("list", length = length(variables))
for(i in seq(from = 0L, to = length(variables))) {
  variables_in_model[[i + 1L]] <- combn(variables, i, simplify = FALSE)
}

variables_in_model <- unlist(variables_in_model, recursive = FALSE)
fit_reduced <- purrr::map(
  variables_in_model,
  ~ fit_linear_regression(biometric, weight, .x)
)

purrr::map_dfr(
  fit_reduced,
  ~ tibble(
    p = length(.x$betas),
    rsq = .x$rsq,
    rsqadj = .x$rsqadj,
    cp = mallows_c(.x, fit),
    mse = .x$mse,
    terms = paste(names(.x$betas), collapse = ", ")
  )
)

위 과정은 함수 evaluate_linear_regression()에 구현되어 있다.

evaluate_linear_regression(biometric, weight, c(age, height))

단계적 변수선택

완전모형에 포함되는 $k$개의 변수 중에서 적절한 독립변수를 단계적으로 선택하는 방법 중 단계별방법(stepwise method)에 대해 알아보자. 이 방법은 $k$개의 독립변수 후보들 중에서 종속변수에 가장 큰 영향을 주는 변수들부터 선택하며 모형에 포함시키면서, 새롭게 모형에 추가된 변수에 기인하여 기존 변수가 그 중요도가 약화되어 모형으로부터 제거될 수 있는지를 매 단계별로 검토하여 해당 변수를 제거하는 과정을 거치며, 추가 또는 제거되는 변수가 더 이상 없을 때 변수 선택을 완료하는 방법이다. 이 때, 각 변수의 추가/제거 단계에서 개별 기울기 회귀계수에 대한 $F$-검정을 적용한다. 이 때, 각 변수를 추가할 때 적용하는 유의수준 $\alpha_{in}$은 제거할 때 적용하는 유의수준 $\alpha_{out}$보다 작거나 같게 설정한다. 즉, 변수가 포함되는 것을 변수가 제거되는 것보다 어렵게 한다.

각각의 변수 추가/제거 단계에서 적용되는 $F$ 값($\alpha_{in}$과 $\alpha_{out}$에 대응하는 값)들을 각각 $F_{in}$, $F_{out}$이라 하자.

alpha_in <- 0.05
alpha_out <- 0.10
variables <- c("age", "height")

• [단계 1] $k$개의 변수 각각을 이용하여 $k$개의 회귀모형을 구하고, $SSR$이 가장 큰 모형에 해당하는 변수($x_{(1)}$이라 하자)를 선택한다. 해당 모형에서 변수에 대한 $F$값이 $F_{in}$보다 크면 변수 $x_{(1)}$을 모형에 포함시킨다. 그렇지 않으면, 아무 변수도 모형에 포함시키지 않은 채 변수 선택을 종료한다.

variables_in_model <- purrr::map_dfr(
  variables,
  ~ test_type2_linear_regression(biometric, weight, .x)
) %>%
  dplyr::slice_max(ss, n = 1L, with_ties = FALSE) %>%
  dplyr::filter(p_value < alpha_in) %>%
  pull(terms)

variables_in_model

• [단계 2] 아직 모형에 포함되지 않은 변수 각각에 대해, 변수를 모형에 추가할 때 해당 변수에 대한 Type $\text{II}$ 제곱합을 산출한다. 최대 제곱합을 갖는 변수($x_j$라 하자)를 선택하고, 그 변수에 해당하는 $F$ 값을 산출한다. $F$ 값이 $F_{in}$보다 크면 변수 $x_j$를 모형에 포함시킨다. 그렇지 않으면, 변수 선택을 종료한다.

variables_not_in_model <- setdiff(variables, variables_in_model)

new_variable_in <- purrr::map_dfr(
  variables_not_in_model,
  ~ test_type2_linear_regression(biometric, weight, c(variables_in_model, .x), 
                                 .last_only = TRUE)
) %>%
  dplyr::slice_max(ss, n = 1L, with_ties = FALSE) %>%
  dplyr::filter(p_value < alpha_in) %>%
  pull(terms)

if (length(new_variable_in) == 0L) {
  print("변수 선택 종료")
} else {
  variables_in_model <- c(variables_in_model, new_variable_in)
  variables_not_in_model <- setdiff(variables_not_in_model, new_variable_in)
}

• [단계 3] 모형에 포함된 변수 각각에 대해 Type $\text{II}$ 제곱합을 산출한다. 최소 제곱합을 갖는 변수($x_i$라 하자)를 선택하고, 그 변수에 해당하는 $F$ 값을 산출한다. $F$ 값이 $F_{out}$보다 작으면 변수 $x_i$를 모형에서 제거하고, 그렇지 않으면 그대로 둔다. [단계 3]의 결과 모든 변수가 모형에 포함되어 있으면 변수 선택을 종료하고, 모형에 포함되지 않은 변수가 남아있는 경우 [단계 2]로 간다.

new_variable_out <- test_type2_linear_regression(biometric, weight, variables_in_model) %>%
  dplyr::slice_min(ss, n = 1L, with_ties = FALSE) %>%
  dplyr::filter(p_value > alpha_out) %>%
  pull(terms)

if (length(new_variable_out) > 0L) {
  variables_in_model <- setdiff(variables_in_model, new_variable_out)
  variables_not_in_model <- c(variables_not_in_model, new_variable_out)
}

if (length(variables_not_in_model) == 0L) {
  print("변수 선택 종료")
}

variables_in_model

위 단계적 변수 선택 방법을 함수 select_variables_stepwise()로 수행할 수 있다.

variables_selected <- select_variables_stepwise(biometric, weight, c(age, height))
variables_selected

위 예제 데이터에 대하여는 모든 변수가 선택된다.

다른 예로, 아래 세 개의 독립변수와 하나의 종속변수로 이루어진 데이터 프레임 pcrdata에 대해 단계별 변수 선택방법에 의해 회귀모형에 사용할 독립변수를 찾아보자.

data(pcrdata, package = "dmtr")
pcrdata

pcrdata_variable_selected <- select_variables_stepwise(pcrdata, y, c(x1, x2, x3))
pcrdata_variable_selected

회귀모형의 진단

잔차산점도

정규확률분포도

반응치에 대한 추정 및 예측

회귀모형이 추정된 후 독립변수의 새로운 관측치에 대하여 종속변수값을 예측한다. 독립변수의 새로운 관측치 $\mathbf{x}{new} = c(1, x{new, 1}, x_{new, 2}, \cdots, x_{new, k})$에 대응하는 종속변수값은 확률변수라 할 수 있는데, 이를 미래반응치(future response) $y_{new}$라 한다. 동일한 독립변수 관측치 $\mathbf{x}{new}$ 수준에서 종속변수값을 여러 번 측정할 때 평균적으로 취하는 값(기대치)를 평균반응치(mean response) $E[y{new} \, | \, \mathbf{x}_{new}]$라 한다.

위 나이/키/몸무게 데이터로 이루어진 회귀모형 예제에서, 나이가 40세이고 키가 170cm이며, 몸무게는 측정이 되지 않은 새로운 관측치가 있다고 하자.

x_new <- c(`(Intercept)` = 1, age = 40, height = 170)

평균반응치의 추정

평균반응치는 아래와 같다.

[ E\left[y_{new} \, | \, \mathbf{x}{new}\right] = \mathbf{x}{new}^{\top} \boldsymbol{\beta} ]

이 때, 평균반응치의 추정량은 위에서 추정한 회귀모형을 이용하여 다음과 같이 구할 수 있다.

[ \hat{y}{new} = \mathbf{x}{new}^{\top} \hat{\boldsymbol{\beta}} ]

y_new_hat <- sum(x_new * fit$betas)
y_new_hat

평균반응치에 대한 신뢰구간을 구하기 위해 우선 $\hat{y}_{new}$의 분산을 아래와 같이 구한다.

[ Var(\hat{y}{new}) = \mathbf{x}{new}^{\top} Var(\hat{\beta}) \mathbf{x}{new} = \sigma^2 \mathbf{x}{new}^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{x}_{new} ]

위 식에 $\sigma^2$ 대신 $MSE$를 대입하여 $\hat{y}_{new}$의 분산을 추정한다.

[ \widehat{Var}(\hat{y}{new}) = MSE \, \mathbf{x}{new}^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{x}_{new} ]

위 분산 추정치에 제곱근을 취하면 $\hat{y}_{new}$의 표준오차를 구할 수 있다.

[ se(\hat{y}{new}) = \sqrt{MSE \, \mathbf{x}{new}^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{x}_{new}} ]

y_new_se <- sqrt(fit$mse * t(x_new) %*% solve(fit$hessian / 2) %*% x_new) %>% drop()
y_new_se

이 때, 평균반응치 $E\left[y_{new} \, | \, \mathbf{x}_{new}\right]$의 $100(1 - \alpha)\%$ 신뢰구간은 다음과 같이 산출된다.

[ \hat{y}{new} \pm t{(\alpha/2, \, n - k - 1)} se(\hat{y}_{new}) ]

alpha <- 0.05
y_new_confint <- c(
  y_new_hat + qt(alpha / 2, df = fit$df) * y_new_se,
  y_new_hat + qt(alpha / 2, df = fit$df, lower.tail = FALSE) * y_new_se
)
y_new_confint

위 과정을 함수 predict_linear_regression()으로 수행할 수 있다. 이 때, .ci_interval은 $(1 - \alpha)$값을 나타낸다.

predict_linear_regression(
  fit,
  .new_data = tibble(age = 40, height = 170),
  .xvar = c(age, height),
  .ci_interval = 0.95
)

미래반응치의 추정

미래반응치 $y_{new}$의 예측치는 평균반응치의 추정치 $\hat{y}_{new}$와 동일하다.

미래반응치의 예측구간을 구하기 위해서는 예측오차의 분산을 알아야 하며, 이는 아래와 같이 수식으로 표현된다.

[ Var(y_{new} - \hat{y}{new}) = \sigma^2 + \sigma^2 \mathbf{x}{new}^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{x}_{new} ]

즉, 예측오차의 분산은 회귀모형의 오차항의 분산 $\sigma^2$과 평균반응치 추정치의 분산 $\sigma^2 \mathbf{x}{new}^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{x}{new}$의 합이다.

이 때, 미래반응치 $y_{new}$의 $100(1 - \alpha)\%$ 예측구간은 다음과 같이 산출된다.

[ \hat{y}{new} \pm t{(\alpha/2, \, n - k - 1)} \sqrt{MSE + se^2(\hat{y}_{new})} ]

alpha <- 0.05
y_new_predint <- c(
  y_new_hat + qt(alpha / 2, df = fit$df) * sqrt(fit$mse + y_new_se ^ 2),
  y_new_hat + qt(alpha / 2, df = fit$df, lower.tail = FALSE) * sqrt(fit$mse + y_new_se ^ 2)
)
y_new_predint

위 과정을 함수 predict_linear_regression()으로 수행할 수 있다. 이 때, .pi_interval은 $(1 - \alpha)$값을 나타낸다.

predict_linear_regression(
  fit,
  .new_data = tibble(age = 40, height = 170),
  .xvar = c(age, height),
  .pi_interval = 0.95
)

예측성능 평가

회귀분석을 예측의 목적으로 사용하는 경우, 예측성능을 평가하는 것이 중요하다. 이를 위해서 학습표본을 훈련표본(training sample)과 테스트 표본(test sample)의 두 부분으로 나누고, 훈련표본만을 사용하여 회귀모형을 구축한 후, 테스트 표본의 독립변수를 사용하여 종속변수를 예측하고, 실제 종속변수값과 비교하여 예측성능을 평가하는 것이다.

예측오차의 척도로 여러 가지를 사용할 수 있으나, 평균절대오차(mean absolute deviation; MAD), 평균제곱오차(mean squared error; MSE) 또는 평균제곱오차의 제곱근(root mean squared error; RMS)을 널리 사용한다.

위에서 회귀모형 추정에 사용한 관측치들 외에 아래와 같이 실제 종속변수값이 관측된 세 개의 테스트 표본이 존재한다고 하자.

test_sample <- tribble(
  ~age, ~height, ~weight,
  30, 175, 75,
  40, 170, 68,
  50, 165, 60
)

우선, 추정된 회귀모형을 이용하여 각 테스트 표본에 대한 예측값을 구하고, 실제 관측값과의 차이를 계산한다.

test_prediction_error <- test_sample[["weight"]] - 
  (predict_linear_regression(fit, test_sample, c(age, height)) %>% pull(.pred))
test_prediction_error

이 때, 평균절대오차, 평균제곱오차 및 그 제곱근은 아래와 같이 계산된다.

tibble(
  mad = mean(abs(test_prediction_error)),
  mse = mean(test_prediction_error ^ 2),
  rmse = sqrt(mean(test_prediction_error ^ 2))
)

위 각각의 평가척도는 함수 eval_mad(), eval_mse(), 그리고 eval_rmse()로 확인할 수 있다. 각 함수들은 관측치 벡터 .y와 예측치 벡터 .y_hat을 입력값으로 사용한다.

list(
  mad = eval_mad,
  mse = eval_mse,
  rmse = eval_rmse
) %>%
  purrr::map_dbl(
    rlang::invoke,
    .args = list(
      .y = test_sample[["weight"]],
      .y_hat = predict_linear_regression(fit, test_sample, c(age, height)) %>% pull(.pred)    
    )
  )

다중공선성

독립변수들 사이에 상당히 높은 선형관계가 존재하는 현상을 다중공선성(multicollinearity)이라고 한다. 완전공선성(perfect collinearity)이 존재하는 경우에는 $\mathbf{X}^{\top} \mathbf{X}$의 역행렬이 존재하지 않아 회귀계수를 추정할 수 없게 된다. 완전공선성은 없다 하여도 다중공선성이 경우에는 추정된 회귀계수의 분산이 매우 커지며, 따라서 정확한 모수추정 및 검정에 어려움이 있다.

다중공선성 척도 - 분산팽창계수(variance inflation factor; $VIF$)

$j$번째 회귀계수의 추정량 $\hat{\beta}_j$에 대한 분산팽창계수 $VIF_j$는 다음과 같이 정의된다.

[ VIF_j = \frac{1}{1 - R_j^2} ]

여기서 $R_j$는 $j$번째 변수를 독립변수로 하고 나머지 $1, \ldots, j -1, j + 1, \ldots, k$ 번째 변수들을 독립변수로 하는 회귀모형에서의 결정계수를 말한다.

각 독립변수에 대한 분산팽창계수를 계산은 함수 vif_linear_regression()으로 구현되어 있다.

vif_linear_regression(biometric, c(age, height))

일반적으로 $k$개의 $VIF_j$ 중 가장 큰 값이 5를 넘으면(가장 큰 $R_j^2$이 0.8을 넘으면) 다중공선성이 있다고 할 수 있으며, 10보다 큰 값이면 심각하다고 볼 수 있다.

다중공선성 해결 방법

다중공선성을 해결하기 위한 이론적 방법으로는 모수추정을 위해 최소자승법 대신에 주성분회귀(principal component regression), 부분최소자승 회귀(partial least squares regression)를 사용하는 방법이 있다. 이 방법들은 다른 장에서 별도로 다룬다.

지시변수와 회귀모형

독립변수가 정량적 변수가 아니라 범주형 변수일 경우, 지시변수를 생성하여 모형을 추정한다. 범주의 수가 $m$인 범주형 변수에 대해서 $m - 1$개의 지시변수를 생성하여야 한다. 이 때, 모형의 분석에 있어서 정량적 변수로만 이루어진 모형과 다르게 고려해야 할 부분들이 있다. 예를 들어, 위에서 살펴본 분산팽창계수의 경우, 범주형 변수가 포함될 때의 분산팽창계수는 보다 일반화된 식을 필요로 한다.

본 패키지에서는 범주형 변수를 이용한 회귀모형은 고려하지 않기로 한다.

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