regression.Rmd
data(biometric, package = "dmtr")| age | height | weight |
|---|---|---|
| 21 | 170 | 60 |
| 47 | 167 | 65 |
| 36 | 173 | 67 |
| 15 | 165 | 54 |
| 54 | 168 | 73 |
| 25 | 177 | 71 |
| 32 | 169 | 68 |
| 18 | 172 | 62 |
| 43 | 171 | 66 |
| 28 | 175 | 68 |
아래와 같이 \(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 \tag{2.1} \end{equation}\]
이 때, 오차항 \(\epsilon_i\)은 서로 독립이고 동일한 정규분포 \(N(0, \sigma^2)\)을 따른다.
위 회귀모형은 아래와 같이 행렬의 연산으로 표한할 수 있다.
\[\begin{equation} \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} \tag{2.2} \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} \]
이다.
다중회귀모형에서 회귀계수들의 추정은 최소자승법(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 \tag{3.1} \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 \tag{3.2} \end{equation}\]
식 (3.2)는 아래와 같이 행렬식으로 표현할 수 있다.
\[ \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)
#> [,1]
#> (Intercept) -108.1671993
#> age 0.3291212
#> height 0.9552913추정된 회귀계수에서 각 편미분값이 0을 확인해보자.
t(X) %*% (y - X %*% beta_hat)
#> [,1]
#> (Intercept) -2.596323e-11
#> age -8.215011e-10
#> height -4.399752e-09실제 얻어진 값은 정확히 0은 아니지만 0에 매우 근접하며, 정확히 0이 얻어지지 않는 이유는 컴퓨팅 측면의 문제로 볼 수 있다 (알고리즘, 소수점 표현 구조 등). 이러한 경우 dplyr::near() 함수를 이용하여 수치가 목표 수치에 충분히 근접하여 있는지를 확인할 수 있으며, 위 회귀계수 편미분값의 경우 모두 0에 충분히 근접함을 확인할 수 있다.
위와 같은 회귀계수 추정과정을 함수 fit_linear_regression()으로 구현하였으며, 결과 리스트값에서 betas 원소가 추정된 회귀계수를 나타낸다.
fit <- fit_linear_regression(biometric, weight, c(age, height))
fit$betas
#> (Intercept) age height
#> -108.1274927 0.3291100 0.9550617위 추정된 회귀계수를 식 (3.1)에 대입하여 잔차제곱합(residual sum of squareds; \(SSE\))을 계산해보자.
sse <- sum((y - (X %*% fit$beta)) ^ 2)
sse
#> [1] 49.26926이 때, 식 (2.1)의 오차항 \(\epsilon_i\)의 분산 \(\sigma ^ 2\)의 추정은 위 잔차제곱합을 자유도 \((n - k - 1)\)로 나눈 잔차평균제곱합(mean sequared error; \(MSE\))을 이용한다.
\[ \hat{\sigma} ^ 2 = MSE = \frac{SSE}{n - k - 1} \]
함수 fit_linear_regression()의 결과 리스트값에서 mse 원소가 잔차평균제곱값을 나타내며, df 원소가 자유도 \((n - k - 1)\)을 나타낸다.
fit$mse
#> [1] 7.038465
fit$df
#> [1] 7따라서, 이 두 결과값을 이용하여 잔차제곱합을 역으로 계산할 수 있다.
fit$mse * fit$df
#> [1] 49.26926다중회귀모형에 있어서 유의성 검정, 즉 회귀성 검정은 독립변수의 기울기에 해당하는 모든 회귀계수(\(\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
#> [1] TRUE위 결과에서 귀무가설은 기각되며, 유의수준 0.05에서 유의한 모형이라 할 수 있다. 위 결과를 함수 anova_linear_regression()을 이용하여 분산분석표로 나타낼 수 있다.
anova_linear_regression(fit)
#> # A tibble: 3 × 6
#> source df ss ms F_statistic p_value
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Model 2 227. 114. 16.1 0.00239
#> 2 Residuals 7 49.3 7.04 NA NA
#> 3 Total 9 276. NA NA NA회귀계수벡터 \(\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)
#> (Intercept) age height
#> (Intercept) 1774.3282944 -0.671107371 -10.264886484
#> age -0.6711074 0.004794718 0.003035477
#> height -10.2648865 0.003035477 0.059566812함수 fit_linear_regression()의 결과 리스트값에서 hessian 원소는 \(2 \mathbf{X}^\top \mathbf{X}\) 값을 지니며, 따라서 hessian 원소를 2로 나누어 \(\mathbf{X}^\top \mathbf{X}\)에 대입할 때 위와 동일한 분산-공분산행렬을 보인다.
fit$mse * solve(fit$hessian / 2)
#> (Intercept) age height
#> (Intercept) 1774.3282947 -0.671107371 -10.264886485
#> age -0.6711074 0.004794718 0.003035477
#> height -10.2648865 0.003035477 0.059566812참고로, hessian 행렬은 식 (3.1)의 오차항 제곱합 \(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)))
#> (Intercept) age height
#> 42.1227764 0.0692439 0.2440631이 표본오차는 fit_linear_regression()의 결과 리스트의 se 원소에 저장되어 있다.
fit$se
#> (Intercept) age height
#> 42.1227764 0.0692439 0.2440631이 때, \(\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
#> (Intercept) age height
#> TRUE TRUE TRUE위 결과에서 모든 각각의 회귀계수에 대해 귀무가설은 기각되며, 유의수준 0.05에서 각 독립변수는 종속변수인 몸무게에 유의한 변수라는 것을 의미한다. 위 결과를 함수 ttest_linear_regression()을 이용하여 데이터 프레임으로 나타낼 수 있다.
ttest_linear_regression(fit)
#> # A tibble: 3 × 5
#> term estimate std_error t_statistic p_value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) -108. 42.1 -2.57 0.0372
#> 2 age 0.329 0.0692 4.75 0.00208
#> 3 height 0.955 0.244 3.91 0.00580\(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
#> [1] 159.0108Type \(\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)
#> [1] TRUE위 결과에서, 유의수준 0.05에서 귀무가설은 기각되며, 따라서 height 변수가 이미 존재하는 모형에 age 변수를 추가할 때 종속변수를 추가로 유의하게 설명한다 할 수 있다.
이와 같은 \(F\) 검정을 회귀모형의 모든 각각의 변수에 수행하는 함수 test_type2_linear_regression()을 아래와 같이 사용할 수 있다.
test_type2_linear_regression(biometric, weight, c(age, height))
#> # A tibble: 2 × 4
#> terms ss F_statistic p_value
#> <chr> <dbl> <dbl> <dbl>
#> 1 age 159. 22.6 0.00208
#> 2 height 108. 15.3 0.00579종속변수를 설명하기 위한 독립변수들의 후보가 매우 많은 경우, 어떤 독립변수들의 조합에 의한 회귀모형이 가장 좋은가를 판단한다.
고려하고 있는 회귀모형 또는 추정된 회귀식이 얼마나 데이터를 잘 반영하고 있는가를 알기 위해, 회귀모형의 결정계수(coefficient of determination) \(R ^ 2\)를 아래와 같이 정의한다.
\[ R ^ 2 = \frac{SSR}{SST} = 1 - \frac{SSE}{SST} \]
즉, \(R ^ 2\)값은 전체제곱합 중 모형이 설명하는 제곱합의 비율로 해석되며, 0에서 1 사이의 값을 갖는다. 위 예제 데이터에서 추정된 회귀모형의 결정계수 값을 계산해보자.
이 식에서, \(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))
#> [1] 0.770817위 결정계수 \(R ^ 2\) 및 수정결정계수 \(R_{adj} ^ 2\)는 아래와 같이 fit_linear_regression() 결과 리스트의 rsq 및 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)
#> [1] 23.59169참고로, \(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)
)
#> Note: Using an external vector in selections is ambiguous.
#> ℹ Use `all_of(.x)` instead of `.x` to silence this message.
#> ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
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 = ", ")
)
)
#> # A tibble: 4 × 6
#> p rsq rsqadj cp mse terms
#> <int> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 1 0 0 31.3 30.7 (Intercept)
#> 2 2 0.432 0.361 16.3 19.6 (Intercept), age
#> 3 2 0.246 0.152 23.6 26.0 (Intercept), height
#> 4 3 0.822 0.771 3 7.04 (Intercept), age, height위 과정은 함수 evaluate_linear_regression()에 구현되어 있다.
evaluate_linear_regression(biometric, weight, c(age, height))
#> # A tibble: 4 × 6
#> p rsq rsqadj cp mse terms
#> <int> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 1 0 0 31.3 30.7 (Intercept)
#> 2 2 0.432 0.361 16.3 19.6 (Intercept), age
#> 3 2 0.246 0.152 23.6 26.0 (Intercept), height
#> 4 3 0.822 0.771 3 7.04 (Intercept), 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")
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
#> [1] "age"
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)
#> Note: Using an external vector in selections is ambiguous.
#> ℹ Use `all_of(variables_in_model)` instead of `variables_in_model` to silence this message.
#> ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
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)
}
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("변수 선택 종료")
}
#> [1] "변수 선택 종료"
variables_in_model
#> [1] "age" "height"위 단계적 변수 선택 방법을 함수 select_variables_stepwise()로 수행할 수 있다.
variables_selected <- select_variables_stepwise(biometric, weight, c(age, height))
#> Inital variable: age , p-value = 0.03902789 < 0.05
#> Iteration 1 : forward selection - height , p-value = 0.00579269 < 0.05
#> Iteration 1 : backward elimination - no variable can be deleted without a statistically significant loss of fit.
#> Every variable provides statistically significant improvement of fit.
variables_selected
#> [1] "age" "height"위 예제 데이터에 대하여는 모든 변수가 선택된다.
다른 예로, 아래 세 개의 독립변수와 하나의 종속변수로 이루어진 데이터 프레임 pcrdata에 대해 단계별 변수 선택방법에 의해 회귀모형에 사용할 독립변수를 찾아보자.
data(pcrdata, package = "dmtr")
pcrdata
#> # A tibble: 6 × 4
#> x1 x2 x3 y
#> <dbl> <dbl> <dbl> <dbl>
#> 1 -3 -3 5 -30
#> 2 -2 -3 7 -20
#> 3 0 0 4 0
#> 4 1 2 0 5
#> 5 2 2 -5 10
#> 6 2 2 -11 35
pcrdata_variable_selected <- select_variables_stepwise(pcrdata, y, c(x1, x2, x3))
#> Inital variable: x1 , p-value = 0.00680743 < 0.05
#> Iteration 1 : forward selection - no additional variable gives statistically significant improvement of the fit.
pcrdata_variable_selected
#> [1] "x1"회귀모형이 추정된 후 독립변수의 새로운 관측치에 대하여 종속변수값을 예측한다. 독립변수의 새로운 관측치 \(\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
#> [1] 67.39739평균반응치에 대한 신뢰구간을 구하기 위해 우선 \(\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
#> [1] 1.006575이 때, 평균반응치 \(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
#> [1] 65.01722 69.77757위 과정을 함수 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
)
#> # A tibble: 1 × 4
#> .pred .se .ci_lower .ci_upper
#> <dbl> <dbl> <dbl> <dbl>
#> 1 67.4 1.01 65.0 69.8미래반응치 \(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
#> [1] 60.68767 74.10712위 과정을 함수 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
)
#> # A tibble: 1 × 4
#> .pred .se .pi_lower .pi_upper
#> <dbl> <dbl> <dbl> <dbl>
#> 1 67.4 1.01 60.7 74.1회귀분석을 예측의 목적으로 사용하는 경우, 예측성능을 평가하는 것이 중요하다. 이를 위해서 학습표본을 훈련표본(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
#> [1] 6.1183987 0.6026067 -5.9131854이 때, 평균절대오차, 평균제곱오차 및 그 제곱근은 아래와 같이 계산된다.
tibble(
mad = mean(abs(test_prediction_error)),
mse = mean(test_prediction_error ^ 2),
rmse = sqrt(mean(test_prediction_error ^ 2))
)
#> # A tibble: 1 × 3
#> mad mse rmse
#> <dbl> <dbl> <dbl>
#> 1 4.21 24.3 4.92위 각각의 평가척도는 함수 eval_mad(), eval_mse(), 그리고 eval_rmse()로 확인할 수 있다. 각 함수들은 관측치 벡터 .y와 예측치 벡터 .y_hat을 입력값으로 사용한다.
독립변수들 사이에 상당히 높은 선형관계가 존재하는 현상을 다중공선성(multicollinearity)이라고 한다. 완전공선성(perfect collinearity)이 존재하는 경우에는 \(\mathbf{X}^{\top} \mathbf{X}\)의 역행렬이 존재하지 않아 회귀계수를 추정할 수 없게 된다. 완전공선성은 없다 하여도 다중공선성이 경우에는 추정된 회귀계수의 분산이 매우 커지며, 따라서 정확한 모수추정 및 검정에 어려움이 있다.
\(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))
#> # A tibble: 2 × 2
#> terms vif
#> <chr> <dbl>
#> 1 age 1.03
#> 2 height 1.03일반적으로 \(k\)개의 \(VIF_j\) 중 가장 큰 값이 5를 넘으면(가장 큰 \(R_j^2\)이 0.8을 넘으면) 다중공선성이 있다고 할 수 있으며, 10보다 큰 값이면 심각하다고 볼 수 있다.