4 장 다중선형회귀모형

반응변수의 변동을 하나의 설명변수만으로 충분하게 설명하는 것은 대부분의 경우 불가능할 것이다. 따라서 반응변수와 관련이 있을 것으로 보이는 여러 개의 설명변수를 모형에 포함시키는 다중회귀모형이 실제 상황에서 많이 사용되는 모형이 된다. 기본 가정 및 추론 방법에서는 단순회귀모형과 큰 차이가 있는 것은 아니지만, 모형에 여러 개의 설명변수가 포함되기 때문에 단순회귀모형에서는 없었던 문제들이 발생할 수 있다.

4.1 다중회귀모형의 설정

다중선형회귀모형에서는 반응변수 \(Y\) 와 설명변수 \(X_{1}, \ldots, X_{k}\) 사이에 다음의 관계가 존재한다고 가정한다.

\[\begin{equation} Y_{i} = \beta_{0} + \beta_{1}X_{1i} + \cdots + \beta_{k}X_{ki} + \varepsilon_{i},~~i=1,\ldots,n \tag{4.1} \end{equation}\]

\(Y_{i}\) 는 반응변수의 \(i\) 번째 값, \(X_{ji}\) 는 설명변수 \(X_{j}\) 의 \(i\) 번째 값을 나타내며, 모수 \(\beta_{0}, \beta_{1}, \ldots, \beta_{k}\) 는 \((k+1)\) 차원의 공간에서 정의되는 회귀평면을 구성하고 있다. 오차항 \(\varepsilon_{i}\) 은 단순회귀모형에서와 동일한 의미를 갖고 있는 확률변수로서, \(\varepsilon_{i} \stackrel{iid}{\sim} N(0,\sigma^{2}),~i=1,\ldots,n\) 의 가정을 하고 있다.

식 (4.1)의 관계식은 다음과 같이 풀어서 표현할 수 있으며,

\[\begin{align*} Y_{1} &= \beta_{0} + \beta_{1}X_{11} + \cdots + \beta_{k}X_{k1} + \varepsilon_{1} \\ Y_{2} &= \beta_{0} + \beta_{1}X_{12} + \cdots + \beta_{k}X_{k2} + \varepsilon_{2} \\ \vdots \\ Y_{n} &= \beta_{0} + \beta_{1}X_{1n} + \cdots + \beta_{k}X_{kn} + \varepsilon_{n} \end{align*}\]

이것은 다시 다음과 같이 행렬의 형태로 정리할 수 있다.

\[ \begin{pmatrix} Y_{1} \\ Y_{2} \\ \vdots \\ Y_{n} \end{pmatrix} = \begin{pmatrix} 1 & X_{11} & \cdots & X_{k1} \\ 1 & X_{12} & \cdots & X_{k2} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & X_{1n} & \cdots & X_{kn} \end{pmatrix} \begin{pmatrix} \beta_{0} \\ \beta_{1} \\ \vdots \\ \beta_{k} \end{pmatrix} + \begin{pmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{n} \end{pmatrix} \tag{4.2} \] 또는 다음의 행렬로 표현할 수 있다.

\[\begin{equation} \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \tag{4.3} \end{equation}\]

\(\mathbf{Y}\) 는 반응변수의 관찰값 벡터이고, \(\mathbf{X}\) 는 설명변수의 관찰값 행렬로써 design matrix라고 하며, \(\boldsymbol{\beta}\) 는 모수 벡터, \(\boldsymbol{\varepsilon}\) 은 오차항 벡터이다.

4.1.1 회귀계수의 추정

회귀계수 \(\beta_{0}, \beta_{1}, \ldots, \beta_{k}\) 의 추정은 단순회귀모형의 경우와 동일한 방식으로 이루어진다. 반응변수 \(Y\) 와 설명변수 \(X_{1}, \ldots, X_{k}\) 에 대해 관측된 \(n\) 개의 자료를 \((y_{i}, x_{1i}, \ldots, x_{ki}), ~i=1,\ldots,n\) 이라고 하면, 다음의 \(RSS\) 를 최소화하는 \(\hat{\boldsymbol{\beta}} = (\hat{\beta}_{0}, \hat{\beta}_{1}, \ldots, \hat{\beta}_{k})\) 를 선택한다.

\[\begin{equation} RSS = \sum_{i=1}^{n} \left(y_{i} - \hat{\beta}_{0} - \hat{\beta}_{1}x_{1i} - \ldots - \hat{\beta}_{k}x_{ki} \right)^{2} \tag{4.4} \end{equation}\]

식 (4.4)의 \(RSS\) 를 최소화시키는 추정값 \(\hat{\beta}_{0}, \ldots, \hat{\beta}_{k}\) 를 구하기 위해서 \(RSS\) 를 \(\hat{\beta}_{j}, ~j = 0, \ldots, k\) 에 대하여 각각 편미분을 실시해서 다음의 \((k+1)\) 개의 방정식을 얻는다.

\[\begin{align} \frac{\partial RSS}{\partial \hat{\beta}_{0}} & = -2 \sum_{i=1}^{n}(y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1}x_{1i}-\cdots-\hat{\beta}_{k}x_{ki}) = 0\\ \frac{\partial RSS}{\partial \hat{\beta}_{1}} & = -2 \sum_{i=1}^{n}x_{1i}(y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1}x_{1i}-\cdots-\hat{\beta}_{k}x_{ki}) = 0 \\ &\vdots \\ \frac{\partial RSS}{\partial \hat{\beta}_{k}} & = -2 \sum_{i=1}^{n}x_{ki}(y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1}x_{1i}-\cdots-\hat{\beta}_{k}x_{ki}) = 0 \tag{4.5} \end{align}\]

식 (4.5)을 정리하면 다음과 같이 표현이 되며,

\[\begin{align} n\hat{\beta}_{0}+\hat{\beta}_{1}\sum_{i=1}^{n} x_{1i} + \cdots + \hat{\beta}_{k}\sum_{i=1}^{n} x_{ki} &= \sum_{i=1}^{n} y_{i} \\ \hat{\beta}_{0}\sum_{i=1}^{n} x_{1i} + \hat{\beta}_{1}\sum_{i=1}^{n} x_{1i}^{2} + \cdots + \hat{\beta}_{k} \sum_{i=1}^{n} x_{1i}x_{ki} &= \sum_{i=1}^{n} x_{1i}y_{i} \\ & \vdots \\ \hat{\beta}_{0}\sum_{i=1}^{n} x_{ki} + \hat{\beta}_{1}\sum_{i=1}^{n} x_{ki}x_{1i} + \cdots + \hat{\beta}_{k}\sum_{i=1}^{n} x_{ki}^{2} &= \sum_{i=1}^{n} x_{ki}y_{i} \tag{4.6} \end{align}\]

행렬의 형태로 표현하면 다음과 같이 된다.

\[ \begin{pmatrix} n & \sum x_{1i} & \cdots & \sum x_{ki} \\ \sum x_{1i} & \sum x_{1i}^{2} & \cdots & \sum x_{1i}x_{ki} \\ \vdots & \vdots & \ddots & \vdots \\ \sum x_{ki} & \sum x_{ki}x_{1i} & \cdots & \sum x_{ki}^{2} \end{pmatrix} \begin{pmatrix} \hat{\beta}_{0} \\ \hat{\beta}_{1} \\ \vdots \\ \hat{\beta}_{k} \end{pmatrix}= \begin{pmatrix} 1 & 1 & \cdots & 1 \\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & \ddots & \vdots \\ x_{k1} & x_{k2} & \cdots & x_{kn} \end{pmatrix} \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{pmatrix} \tag{4.7} \]

식 (4.7)의 좌변 첫 번째 행렬은 식 (4.3)에서 정의된 행렬 \(\mathbf{X}\) 로 다음과 같이 표현된다.

\[ \begin{pmatrix} n & \sum x_{1i} & \cdots & \sum x_{ki} \\ \sum x_{1i} & \sum x_{1i}^{2} & \cdots & \sum x_{1i}x_{ki} \\ \vdots & \vdots & \ddots & \vdots \\ \sum x_{ki} & \sum x_{ki}x_{1i} & \cdots & \sum x_{ki}^{2} \end{pmatrix}= \mathbf{X}^{T}\mathbf{X} \tag{4.8} \] 단, \(\mathbf{X}^{T}\) 는 행렬 \(\mathbf{X}\) 의 전치행렬이다.

따라서 다중회귀모형의 회귀계수에 대한 최소제곱추정량을 구하기 위한 방정식은 다음과 같이 주어진다.

\[\begin{equation} \mathbf{X}^{T}\mathbf{X}\hat{\boldsymbol{\beta}} = \mathbf{X}^{T}\mathbf{Y} \tag{4.9} \end{equation}\]

만일 \(\mathbf{X}^{T}\mathbf{X}\) 의 역행렬이 존재한다면, 최소제곱추정량은 다음과 같이 표현된다.

\[\begin{equation} \hat{\boldsymbol{\beta}} = (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{Y} \tag{4.10} \end{equation}\]

\(\bullet\) 예제: 다중회귀모형의 회귀계수 추정

데이터 프레임 mtcars는 1974년 \(\textit{Motor Trend US}\) 에 실린 32 종류 자동차 모델의 연비와 관련된 자료가 입력되어 있다.

str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

연비를 나타내는 mpg를 반응변수로 하고, 설명변수에는 차량의 무게를 나타내는 wt와 0.25 마일까지 도달하는 데 걸리는 시간을 나타내는 qsec을 선택해서 다중회귀모형을 설정하고 회귀계수를 추정해 보자.

단순회귀모형의 경우와 동일하게 다중회귀모형에서도 함수 lm()을 사용해서 모형적합을 실시할 수 있다. 하지만 단순회귀모형의 경우와 다르게 설명변수의 개수가 2개 이상이 되며, 복잡한 형태의 모형이 설정되는 경우도 많이 있기 때문에 다중회귀모형의 경우에는 함수 lm()에 지정되는 R 모형공식에 반응변수와 설명변수를 구분하는 ~ 기호 외에도 많은 기호가 사용된다.

R 모형공식에서 사용되는 기호
기호	사용법
물결표 (`~`)	반응변수와 설명변수의 구분. 물결표의 왼쪽에는 반응변수, 오른쪽에는 설명변수를 둔다.
플러스 (`+`)	모형에 포함된 설명변수의 구분. 반응변수 `y`와 설명변수 `x1`, `x2`, `x3`의 회귀모형은 `y ~ x1 + x2 + x3`로 표현된다.
콜론 (`:`)	설명변수 사이의 상호작용 표현. 반응변수 `y`와 설명변수 `x1`, `x2` 그리고 `x1`과 `x2`의 상호작용이 포함된 모형은 `y ~ x1 + x2 + x1:x2`로 표현된다.
별표 (`*`)	주효과와 상호작용 효과를 포함한 모든 효과 표현. `y ~ x1 * x2 * x3`는 `y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3 + x1:x2:x3`를 의미한다.
윗꺽쇠 (`^`)	지정된 차수까지의 상호작용 표현. `y ~ (x1 + x2 + x3)^2`는 `y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3`을 의미한다.
점 (`.`)	반응변수를 제외한 데이터 프레임에 있는 모든 변수. 만일 데이터 프레임에 `y`, `x1`, `x2`, `x3`가 있다면 `y ~ .` 은 `y ~ x1 + x2 + x3`을 의미한다.
마이너스(`-`)	회귀모형에서 제외되는 변수. `y ~ (x1 + x2 + x3)^2 – x2:x3`는 `y ~ x1 + x2 + x3 + x1:x2 + x1:x3`을 의미한다.
`- 1`, `+ 0`	절편 제거. `y ~ x – 1` 혹은 `y ~ x + 0`은 원점을 지나는 회귀모형을 의미한다.
`I()`	괄호 안의 연산자를 수학 연산자로 인식. `y ~ x1 + I(x2+x3)`는 \(y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}(x_{2}+x_{3})+\varepsilon\) 모형을 의미한다.

이제 mpg를 반응변수로, wt와 qsec을 설명변수로 하는 다중회귀모형을 적합해 보자.

fit_m <- lm(mpg ~ wt + qsec, data = mtcars)
fit_m
## 
## Call:
## lm(formula = mpg ~ wt + qsec, data = mtcars)
## 
## Coefficients:
## (Intercept)           wt         qsec  
##     19.7462      -5.0480       0.9292

적합된 모형식은 \(\widehat{\mbox{mpg}} = 19.74 - 5.047\mbox{ wt} + 0.929\mbox{ qsec}\) 임을 알 수 있다. 추정된 회귀계수에 대한 해석은 단순회귀모형의 경우와 비슷하지만 제한 사항이 추가된다. 단순회귀모형에서는 하나의 설명변수만 있기 때문에 추정된 모형은 직선으로 표현되고, 따라서 해당 설명변수가 한 단위 증가했을 때 반응변수의 평균 변화량으로 회귀계수를 해석할 수 있었다. 하지만 다중회귀모형에서는 추정된 회귀모형이 직선이 아닌 2차원 이상에서 정의되는 평면이 되기 때문에 조금 더 복잡한 상황이 되는데, 그것은 어느 한 설명변수의 효과를 측정하기 위해서는 회귀평면을 구성하고 있는 다른 설명변수의 값이 고정되어야 하기 때문이다.

그림 4.1에서는 mpg를 반응변수로, wt와 qsec을 설명변수로 하는 다중회귀모형으로 적합된 회귀평면이 작성되어 있다. 변수 wt의 추정된 회귀계수 \(-5.047\) 은 qsec을 일정한 수준으로 고정한 상태에서 wt가 한 단위 증가했을 때 mpg의 평균 변화량을 나타내는 것인데, 그림에서 회귀평면이 wt 축의 양의 방향으로 급격하게 기우는 모습을 볼 수 있다. 또한 변수 qsec에 대한 추정된 회귀계수 \(0.929\) 은 wt를 일정한 수준으로 고정한 상태에서 qsec이 한 단위 증가했을 때 mpg의 평균 변화량을 나타내는 것이며, 그림에서 회귀평면이 qsec의 양의 방향으로 완만하게 증가하는 모습을 볼 수 있다.

그림 4.1: 설명변수가 2개인 다중회귀모형에서 추정된 회귀평면

\(\bullet\) 예제: 행렬 state.x77

행렬 state.x77은 미국 50개 주와 관련된 8개 변수로 구성되었다. 변수 Life Exp와 HS Grad는 이름 중간에 빈칸이 있으며, 각 주의 이름이 행렬의 row name으로 입력되어 있음을 알 수 있다.

state.x77[1:5,]
##            Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
## Alabama          3615   3624        2.1    69.05   15.1    41.3    20  50708
## Alaska            365   6315        1.5    69.31   11.3    66.7   152 566432
## Arizona          2212   4530        1.8    70.55    7.8    58.1    15 113417
## Arkansas         2110   3378        1.9    70.66   10.1    39.9    65  51945
## California      21198   5114        1.1    71.71   10.3    62.6    20 156361

변수 Murder를 반응변수로, 나머지 7개 변수를 설명변수로 하는 다중회귀모형을 설정해 보자. 회귀분석을 원활하게 수행하기 위해서는 자료가 행렬이 아닌 데이터 프레임의 형태로 입력되어 있어야 한다. 행렬 state.x77을 데이터 프레임으로 변환하면서, 빈칸이 있는 변수 이름을 수정하고 반응변수를 마지막 변수로 이동시키자. 반응변수의 위치를 마지막으로 이동시킨 이유는 산점도 행렬을 작성하고 결과를 해석할 때 더 편리하기 때문이다.

states <- as.data.frame(state.x77) |>  
  rename(Life_Exp = `Life Exp`, HS_Grad = `HS Grad`) |>  
  relocate(Murder, .after = last_col())

str(states)
## 'data.frame':    50 obs. of  8 variables:
##  $ Population: num  3615 365 2212 2110 21198 ...
##  $ Income    : num  3624 6315 4530 3378 5114 ...
##  $ Illiteracy: num  2.1 1.5 1.8 1.9 1.1 0.7 1.1 0.9 1.3 2 ...
##  $ Life_Exp  : num  69 69.3 70.5 70.7 71.7 ...
##  $ HS_Grad   : num  41.3 66.7 58.1 39.9 62.6 63.9 56 54.6 52.6 40.6 ...
##  $ Frost     : num  20 152 15 65 20 166 139 103 11 60 ...
##  $ Area      : num  50708 566432 113417 51945 156361 ...
##  $ Murder    : num  15.1 11.3 7.8 10.1 10.3 6.8 3.1 6.2 10.7 13.9 ...

다중회귀모형에 의한 분석을 진행하기 전에 반드시 거쳐야 할 단계가 있는데, 그것은 모형에 포함될 변수들의 관계를 탐색하는 것이다. 변수들 사이의 관계 탐색에 많이 사용되는 방법은 상관계수와 산점도행렬이다. 산점도행렬과 같은 그래프에 의한 탐색은 필수적이라 할 수 있으며, 연속형 변수의 경우에는 두 변수끼리 짝을 지어 상관계수를 살펴보는 것도 필요하다.

상관계수는 변수들 사이의 선형관계 정도를 확인할 때 사용할 수 있는데, 함수 cor()로 계산할 수 있다. 사용법은 다음과 같다.

cor(x, y = NULL, use = "everything", method = c("pearson", "kendall", "spearman"))

x와 y 모두 벡터, 행렬 혹은 데이터 프레임이 가능한데, x만 주어지면 x에 포함된 모든 변수들 사이의 상관계수를 계산하게 되고, y도 주어지면 x에 속한 변수와 y에 속한 변수들을 하나씩 짝을 지어 상관계수를 계산하게 된다.

use는 결측값의 처리방식에 대한 것으로 선택할 수 있는 것은 디폴트인 "everything"과 "all.obs", "complete.obs", "pairwise.complete.obs"이 있으며, 해당 문자의 약칭으로도 사용이 가능하다. 결측값이 존재하면 use = "everything"인 경우에는 NA가 계산결과로 출력되고, use = "all"인 경우에는 오류가 발생한다. 또한 use = "complete"의 경우에는 NA가 있는 행은 모두 제거된 상태에서 상관계수가 계산되고, use = "pairwise"의 경우에는 상관계수가 계산되는 변수들만을 대상으로 NA가 있는 행을 제거하고 상관계수를 계산한다.

method는 계산하는 상관계수의 종류를 선택한다. 디폴트인 "pearson"은 Pearson의 상관계수를 지정하는 것으로 두 변수 사이의 선형관계 정도를 표현하는 가장 일반적으로 많이 사용되는 상관계수를 계산한다. 두 번째 방법인 "kendall"은 Kendall의 순위상관계수 혹은 Kendall의 \(\tau\) 를 지정하는 것으로서 concordant pair와 discordant pair를 이용하여 정의되는 비모수 상관계수를 계산하며, 순서형 자료에 주로 적용되는 방법이다. 세 번째 방법인 "spearman"은 Spearman의 순위상관계수 혹은 Spearman의 \(\rho\) 를 지정하는 것으로서 두 변수 사이의 관계가 단조증가 혹은 단조감소 함수로 얼마나 잘 설명될 수 있는지를 표현하는 비모수 상관계수를 계산한다. 정규성 가정이 어긋나는 경우에 사용할 수 있다.

데이터 프레임 states를 구성하고 있는 변수들의 상관계수를 구해보자.

cor(states)
##             Population     Income  Illiteracy    Life_Exp     HS_Grad
## Population  1.00000000  0.2082276  0.10762237 -0.06805195 -0.09848975
## Income      0.20822756  1.0000000 -0.43707519  0.34025534  0.61993232
## Illiteracy  0.10762237 -0.4370752  1.00000000 -0.58847793 -0.65718861
## Life_Exp   -0.06805195  0.3402553 -0.58847793  1.00000000  0.58221620
## HS_Grad    -0.09848975  0.6199323 -0.65718861  0.58221620  1.00000000
## Frost      -0.33215245  0.2262822 -0.67194697  0.26206801  0.36677970
## Area        0.02254384  0.3633154  0.07726113 -0.10733194  0.33354187
## Murder      0.34364275 -0.2300776  0.70297520 -0.78084575 -0.48797102
##                 Frost        Area     Murder
## Population -0.3321525  0.02254384  0.3436428
## Income      0.2262822  0.36331544 -0.2300776
## Illiteracy -0.6719470  0.07726113  0.7029752
## Life_Exp    0.2620680 -0.10733194 -0.7808458
## HS_Grad     0.3667797  0.33354187 -0.4879710
## Frost       1.0000000  0.05922910 -0.5388834
## Area        0.0592291  1.00000000  0.2283902
## Murder     -0.5388834  0.22839021  1.0000000

함수 cor()에 데이터 프레임을 입력하면 모든 변수들 사이의 상관계수가 행렬 형태로 계산되어 출력된다. 상관계수행렬은 변수의 개수가 많아지면 변수 사이의 관계 파악이 어려워지는 문제가 있다. 이런 경우에는 그래프로 표현하는 것이 더 효과적인 방법이 될 것이다.

패키지 GGally의 함수 ggcorr()은 상관계수행렬을 그래프로 표현할 때 많이 사용되는 함수이며, 사용법은 다음과 같다.

ggcorr(data, method = c("pairwise", "pearson"), label = FALSE, label_round = 1, ...)

method는 결측값 처리 방식과 계산되는 상관계수의 종류를 차례로 지정하는데, 디폴트는 "pairwise"와 "pearson"이다. label은 그래프에 상관계수를 표시할 것인지 여부를 결정하는 것이고, label_round는 표시되는 상관계수의 반올림 자릿수를 지정한다.

데이터 프레임 states를 구성하고 하는 변수들의 상관계수를 그래프로 나타내 보자. 숫자로만 구성되어 있는 상관계수행렬보다 훨씬 간편하게 변수 사이의 상관관계를 파악할 수 있다.

library(GGally)
ggcorr(states, label = TRUE, label_round = 2)

그림 4.2: 함수 ggcorr()에 의한 상관계수 그래프

반응변수를 마지막 변수로 이동시켰기 때문에 반응변수와 다른 설명변수 사이의 상관계수를 편하게 확인할 수 있다. 반응변수 Murder가 설명변수 Illiteracy와는 비교적 높은 양의 상관관계를, Life_Exp와는 비교적 높은 음의 상관관계를 보이고 있다.

상관계수는 두 변수 사이의 선형관계 정도만을 측정하는 측도이다. 변수 사이에 존재하는 ’있는 그대로’의 관계를 확인하는 가장 좋은 방법은 산점도행렬이다. 함수 GGally::ggpairs()로 작성해 보자.

ggpairs(states, lower = list(continuous = "smooth"))

그림 4.3: 함수 ggpairs()에 의한 산점도행렬

반응변수를 마지막 변수로 위치를 이동시켰기 때문에 산점도행렬의 마지막 행을 이루고 있는 모든 패널에서 Murder가 Y축 변수로 배치되고, 설명변수들이 각각 X축 변수로 배치되었다.

이제 states 자료에 대한 회귀모형을 함수 lm()으로 적합해 보자.

fit_s <- lm(Murder ~ ., data = states)

fit_s
## 
## Call:
## lm(formula = Murder ~ ., data = states)
## 
## Coefficients:
## (Intercept)   Population       Income   Illiteracy     Life_Exp      HS_Grad  
##   1.222e+02    1.880e-04   -1.592e-04    1.373e+00   -1.655e+00    3.234e-02  
##       Frost         Area  
##  -1.288e-02    5.967e-06

함수 lm()으로 생성된 객체 fit_s을 단순하게 출력시키면 추정된 회귀계수만 나타나지만, 사실 객체 fit_s는 많은 양의 정보가 담겨 있는 리스트이며, 이것은 획득하기 위해서는 몇 가지 함수를 사용해야 한다. 함수 목록은 아래 표에 정리되어 있으며, 사용하는 방법은 앞으로 살펴보겠다.

함수 `lm()` 으로 생성된 객체에서 필요한 결과를 얻기 위해 유용하게 사용되는 함수
함수	산출 결과
`anova()`	추정된 회귀모형의 분산분석표 혹은 두 개 이상의 추정된 모형을 비교하기 위한 분산분석표
`coefficients()`	추정된 회귀계수. `coef()`도 가능.
`confint()`	회귀계수의 신뢰구간. 95% 신뢰구간이 디폴트.
`deviance()`	잔차제곱합(residual sum of squares; RSS), \(\sum (y_{i}-\hat{y}_{i})^{2}\)
`fitted()`	반응변수의 적합값, \(\hat{y}_{i}\)
`residuals()`	회귀모형의 잔차, \(e_{i}=y_{i}-\hat{y}_{i}\) . `resid()`도 가능.
`summary()`	회귀모형의 다양한 적합 결과

4.1.2 다항회귀모형

식 (4.1)에 정의된 다중회귀모형에서는 반응변수와 설명변수 사이의 관계가 선형이라고 가정하고 있다. 하지만 두 변수의 관계가 그림 1.4와 같이 명확한 2차 함수 관계가 있을 경우에는 설명변수의 제곱항을 모형에 추가하는 다항회귀모형이 더 적절한 대안이 될 수 있다.

단순회귀모형에 대한 \(p\) 차 다항회귀모형은 다음과 같이 설정된다.

\[\begin{equation} Y = \beta_{0} + \beta_{1}X + \beta_{2}X^{2} + \cdots + \beta_{p}X^{p} + \varepsilon \tag{4.11} \end{equation}\]

차수 \(p\) 는 가능한 낮게 잡는 것이 좋은데, 그것은 너무 높은 차수의 항을 모형에 포함시키면 \(\left( \mathbf{X}^{T} \mathbf{X} \right)^{-1}\) 가 불안정해질 수 있고, 따라서 회귀계수 추정에 문제가 생길 수 있기 때문이다. 이 문제는 ’다중공선성’을 소개할 때 다시 살펴보겠다.

일단 차수가 선택되면, 선택된 차수 이하의 모든 차수는 반드시 모형에 포함되어야 한다. 예를 들어 \(Y=\beta_{0}+\beta_{1}X+\beta_{2}X^{2}+\varepsilon\) 모형에서 어떠한 이유에서 \(X\) 의 1차항을 제거한다면, 모형은 \(Y=\beta_{0}+\beta_{2}X^{2}+\varepsilon\) 가 되는데, 이 모형은 Y축을 중심으로 좌우대칭을 이루어야 하는 지나치게 강한 가정을 갖게 된다. 또한 만일 변수 \(X\) 의 값을 \(a\) 만큼 평행 이동시켜야 한다면, \(X\) 가 \(X+a\) 로 변경되어 다음과 같이 모형에 다시 1차항이 나타나게 된다.

\[\begin{align*} Y &= \beta_{0} + \beta_{2}(X+a)^{2} + \varepsilon \\ &= \beta_{0} + \beta_{2}(X^{2} + 2aX + a^{2}) + \varepsilon \end{align*}\]

설명변수가 \(X_{1}\) 과 \(X_{2}\) 인 다중회귀모형에서 두 변수가 모두 반응변수와 2차 함수 관계가 있다면, 다음과 같은 다항회귀모형이 적절한 회귀모형이 될 수 있다.

\[\begin{equation} Y = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \beta_{3}X_{1}^{2} + \beta_{4}X_{2}^{2} + \varepsilon \tag{4.12} \end{equation}\]

\(\bullet\) 예제: women

데이터 프레임 women에는 미국 30대 여성의 키(height)와 몸무게(weight)가 입력되어 있다. 두 변수 height와 weight의 회귀모형을 설정해 보자. 먼저 두 변수의 관계가 선형인지 여부를 그래프로 확인해 보자.

ggplot(women, aes(x = height, y = weight)) + 
   geom_point(size = 3) +
   geom_smooth(aes(color = "loess"), se = FALSE) +
   geom_smooth(aes(color = "linear"), method="lm", se = FALSE) +
   labs(color = NULL)

그림 4.4: 데이터 프레임 women의 변수 height와 weight의 산점도와 선형회귀직선 및 국소회귀곡선

두 변수의 관계가 선형보다는 2차 곡선이 더 적합한 것으로 보인다.

다항회귀모형은 함수 poly()를 이용하거나 또는 함수 I()를 이용해서 적합할 수 있다. 함수 poly()는 함수 lm()과 함께 lm(y ~ poly(x, degree = 1, raw = FALSE), data)와 같이 사용할 수 있다. degree는 다항회귀모형의 차수를 지정하는 것으로 디폴트 값은 1차이고, raw는 직교다항회귀(orthogonal polynomial regression)의 사용 여부를 선택하는 것으로서 디폴트 값인 FALSE는 직교다항회귀에 의한 적합이 된다. 따라서 일반적인 다항회귀모형을 사용하고자 한다면 반드시 raw에 TRUE를 지정해야 한다.

차수가 높지 않은 일반적인 다항회귀모형을 적합하는 경우에는 함수 I()를 사용하는 것이 더 간편할 수 있다. 함수 lm() 안에서 lm(y ~ x + I(x^2), data)와 같이 입력하면 2차 다항회귀모형을 적합하게 된다.

반응변수 weight에 대한 설명변수 height의 2차 다항회귀모형을 적합해 보자.

fit_w <- lm(weight ~ height + I(height^2), women)

fit_w
## 
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
## 
## Coefficients:
## (Intercept)       height  I(height^2)  
##   261.87818     -7.34832      0.08306

적합된 회귀모형식은 \(\hat{y}_{i}=261.87 - 7.348x_{i} + 0.08x_{i}^{2}\) 임을 알 수 있다.

4.1.3 가변수 회귀모형

선형회귀모형에서 반응변수는 유형이 반드시 연속형이어야 하며, 정규분포의 가정이 필요하다. 반면에 설명변수는 연속형 변수와 범주형 변수가 모두 가능한데, 연속형인 경우에는 정규분포의 가정은 필요 없으나 가능한 좌우대칭에 가까운 분포를 하는 것이 좋다. 범주형 변수 중 순서형 변수인 경우에는 연속형 변수처럼 변수의 값을 그대로 사용하는 것이 가능하지만, 명목형 변수의 경우에는 변수의 값을 그대로 사용하면 절대로 안 되고, 반드시 가변수(dummy variable)를 대신 사용해야 한다.

가변수를 설정하는 방식은 몇 가지가 있는데, 그 중 회귀계수 \(\beta_{0}\) 를 유지하는 방식을 살펴보자. 이 경우, 가변수는 0 또는 1의 값을 갖게 되며, 범주의 개수보다 하나 작은 개수의 가변수를 사용하게 된다. 예를 들어, “Yes”, “No”와 같은 2개의 범주를 갖는 범주형 변수와 연속형 변수 \(X\) 를 설명변수로 하는 회귀모형은 다음과 같이 하나의 가변수를 갖게 된다.

\[\begin{equation} Y = \beta_{0} + \beta_{1}X + \beta_{2}D + \varepsilon \tag{4.13} \end{equation}\]

가변수 \(D\) 가 “Yes” 범주이면 1, “No” 범주이면 0을 값으로 갖는다면, 모집단 회귀직선은 다음과 같이 주어진다. \[ E(Y) = \begin{cases} \beta_{0} + \beta_{1}X & \quad \text{if } ~~D = 0 \\ \beta_{0} + \beta_{2} + \beta_{1}X & \quad \text{if } ~~D = 1 \end{cases} \]

\(D=0\) 이 되는 범주를 ’기준범주’라고 하는데, 절편 \(\beta_{0}\) 가 기준범주의 효과를 나타내고 있고, 가변수의 회귀계수 \(\beta_{2}\) 는 기준범주와 “Yes” 범주의 효과 차이를 나타내고 있다.

식 (4.13)의 회귀모형에서는 반응변수 \(Y\) 와 설병변수 \(X\) 가 “Yes”와 “No” 범주에서 동일한 관계를 갖고 있다고 가정하고 있는데, 만일 각 범주에서 두 변수의 관계가 다를 수 있다면 기울기도 다르게 설정할 필요가 있다. 만일 기울기가 그룹별로 다르다고 가정한다면 가변수와 다른 설명변수의 상호작용 효과를 추가해야 한다. 식 (4.13)의 회귀모형에 가변수와 연속형 변수의 상호작용 효과가 추가된 모형은 다음과 같이 표현된다.

\[\begin{equation} Y = \beta_{0} + \beta_{1}X + \beta_{2}D + \beta_{3}DX + \varepsilon \tag{4.14} \end{equation}\]

이 경우에 모집단 회귀직선은 다음과 같다.

\[ E(Y) = \begin{cases} \beta_{0} + \beta_{1}X & \quad \text{if } ~~D = 0 \\ (\beta_{0} + \beta_{2}) + (\beta_{1}+\beta_{3})X & \quad \text{if } ~~D = 1 \end{cases} \]

연속형 설명변수와 가변수의 상호작용 효과를 모형에 추가하는 것이 더 포괄적인 모형을 구축하는 방법이지만, 설명변수가 많은 경우에는 지나치게 많은 변수가 모형에 추가될 수 있어서 분석에 큰 부담이 될 수 있다. 반드시 필요한 상호작용 효과를 선별해서 추가하는 것이 바람직하다.

범주형 변수가 3개의 범주를 갖는 경우에는 2개의 가변수가 필요하게 된다. 예를 들어 “high”, “medium”, “low”의 3가지 범주를 갖는 범주형 변수와 연속형 변수 \(X\) 를 설명변수로 하는 회귀모형은 다음과 같은 형태를 갖게 된다. 반응변수 \(Y\) 와 설병변수 \(X\) 가 세 범주에 동일한 관계를 갖고 있다고 가정하자.

\[\begin{equation} Y = \beta_{0} + \beta_{1}X + \beta_{2}D_{1} + \beta_{3}D_{2} + \varepsilon \end{equation}\]

단, \(D_{1}\) 은 범주가 “high”이면 1, 아니면 0이고, \(D_{2}\) 는 범주가 “medium”이면 1, 아니면 0이 되는 가변수이다. 이 경우, 기준범주는 “low”가 되며, \(\beta_{2}\) 는 기준범주와 “high” 범주의 효과 차이를, \(\beta_{3}\) 은 기준범주와 “medium” 범주의 효과 차이를 나타내고 있다.

\(\bullet\) 예제: mtcars

반응변수는 자동차 연비인 mpg로 하고, 설명변수로는 1/4 마일 도착 시간인 qsec, 변속기 종류를 나타내는 am과 실린더 개수를 나타내는 cyl을 사용하자.

변수 am과 cyl은 범주형 변수이지만 숫자로 입력된 상태이다.

str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

함수 lm()은 요인으로 설명변수를 입력하면 자동으로 필요한 개수의 가변수를 생성한다. 따라서 변수 am과 cyl은 요인으로 유형을 변경하는 것이 필요하다.

mtcars_3 <- mtcars |> 
  mutate(am = factor(am), cyl = factor(cyl))

요인으로 변경된 두 변수의 범주별 도수를 확인해 보자.

mtcars_3 |> count(am)
##   am  n
## 1  0 19
## 2  1 13

mtcars_3 |> count(cyl)
##   cyl  n
## 1   4 11
## 2   6  7
## 3   8 14

두 개의 범주를 갖고 있는 am과 숫자형 변수 qsec을 설명변수로 하는 회귀모형을 적합해 보자.

fit_m1 <- lm(mpg ~ qsec + am, data = mtcars_3)

적합 결과는 다음과 같다.

fit_m1
## 
## Call:
## lm(formula = mpg ~ qsec + am, data = mtcars_3)
## 
## Coefficients:
## (Intercept)         qsec          am1  
##     -18.889        1.982        8.876

적합된 회귀모형을 그래프로 표현하면 그림 4.5과 같이 두 개의 평행한 회귀직선으로 표현된다.

그림 4.5: 모형 fit_m1의 적합 결과

변수 am의 두 범주 중 기준 범주는 am = 0인 범주이며, 모형 fit_m1에서 사용된 가변수 am1은 두 범주의 효과 차이를 나타내고 있는데, 회귀계수가 8.876으로 계산된 것은 am = 0 범주보다 am = 1 범주에서 반응변수의 평균이 8.876 높다고 추정된 것이다.

또한 적합된 두 직선이 평행한 것은 변수 qsec와 mpg의 관계가 am = 0 범주와 am = 1 범주에서 동일하다고 가정했기 때문이다. 동일한 관계를 유지한다는 가정을 하지 않은 상태에서 회귀모형을 적합한 결과는 그림 4.6와 같다. 어떤 모형을 사용하는 것이 더 적절하다고 생각하는가? 이 문제는 회귀모형의 추론을 통해 결정해야 하는 문제가 된다.

그림 4.6에 표시된 회귀모형에는 변수 qsec과 am의 상호작용 항이 포함되었으며, 자세한 설명은 4.2절에서 살펴보겠다.

그림 4.6: 변수 qsec과 am의 상호작용 항이 포함된 회귀모형의 적합 결과

이번에는 세 개의 범주를 갖고 있는 cyl과 숫자형 변수 qsec을 설명변수로 하는 회귀모형을 적합해 보자.

fit_m2 <- lm(mpg ~ qsec + cyl, data = mtcars_3)

fit_m2
## 
## Call:
## lm(formula = mpg ~ qsec + cyl, data = mtcars_3)
## 
## Coefficients:
## (Intercept)         qsec         cyl6         cyl8  
##     35.0724      -0.4394      -7.4305     -12.6029

적합된 모형은 다음과 같이 그래프로 표현된다.

그림 4.7: 모형 fit_m2의 적합 결과

회귀모형 fit_m2에서 cyl의 기준범주는 cyl = 4인 범주이며, 회귀계수의 추정 결과로 기준범주에서 반응변수의 평균이 cyl = 6인 범주보다 7.43 높고, cyl = 8인 범주보다 12.6 높음을 알 수 있다.

4.2 다중회귀모형의 추론

4.1절에서 우리는 다항회귀모형과 가변수 회귀모형 등 몇 가지 형태의 다중회귀모형을 살펴보았고, 최소제곱추정량에 의한 회귀계수의 추정 방법도 살펴보았다. 지금부터는 회귀계수의 추론에 대해 살펴보도록 하자. 반응변수와 설명변수에 대해 자료가 관측되면 회귀모형에 대한 다양한 추론을 할 수 있는데, 이러한 추론을 진행하기 위해서는 회귀계수의 추정량에 대한 표본분포를 알아야 한다.

4.2.1 회귀계수 추정량의 표본분포

다중회귀모형은 다음과 같이 행렬의 형태로 나타내는 것이 일반적인 표현 방식이 된다.

\[\begin{equation} \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \end{equation}\]

반응변수 벡터 \(\mathbf{Y}\) 와 설명변수 관찰값 행렬 \(\mathbf{X}\), 모수 벡터 \(\boldsymbol{\beta}\), 오차항 벡터 \(\boldsymbol{\varepsilon}\) 의 정의는 식 (4.2) 에서 볼 수 있다. 또한 오차항 \(\varepsilon_{i}, ~i=1,\ldots,n\) 이 모두 평균이 \(0\), 분산이 \(\sigma^{2}\) 이며, 서로 독립이라는 가정도 다음과 같이 행렬로 표현할 수 있다.

\[ E(\boldsymbol{\varepsilon}) = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \]

\[ Var(\boldsymbol{\varepsilon}) = \begin{pmatrix} \sigma^{2} & 0 & \cdots & 0 \\ 0 & \sigma^{2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma^{2} \end{pmatrix}= \sigma^{2}~\mathbf{I} \]

단, \(\mathbf{I}\) 는 \(n\) 차 단위행렬이다.

오차항 벡터의 가정을 기반으로 반응변수 벡터 \(\mathbf{Y}\) 의 평균과 분산은 다음과 같이 구할 수 있으며,

\[\begin{equation} E(\mathbf{Y}) = E(\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}) = E(\mathbf{X}\boldsymbol{\beta}) = \mathbf{X}\boldsymbol{\beta} \end{equation}\]

\[\begin{equation} Var(\mathbf{Y}) = Var(\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}) = Var(\boldsymbol{\varepsilon}) = \sigma^{2}~\mathbf{I} \end{equation}\]

오차항이 정규분포한다는 가정을 추가하면 반응변수 벡터 \(\mathbf{Y}\) 의 분포는 다음과 같다.

\[\begin{equation} \mathbf{Y} \sim N(\boldsymbol{\beta},~\sigma^{2}\mathbf{I}) \end{equation}\]

이제 모수 벡터 \(\boldsymbol{\beta}\) 의 최소제곱추정량 \(\hat{\boldsymbol{\beta}} = (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{Y}\) 의 표본분포를 알아보자. 먼저 추정량 벡터의 평균은 다음과 같이 구할 수 있다.

\[\begin{align*} E(\hat{\boldsymbol{\beta}}) &= E((\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{Y}) \\ &= (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}E(\mathbf{Y}) \\ &= (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{X}\boldsymbol{\beta} \\ &= \boldsymbol{\beta} \end{align*}\]

추정량 벡터의 분산은 다음과 같이 구할 수 있다.

\[\begin{align*} Var(\hat{\boldsymbol{\beta}}) &= Var((\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{Y}) \\ &= (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}~Var(\mathbf{Y})~ ((\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T})^{T} \\ &= (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}~\sigma^{2}\mathbf{I}~ ((\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T})^{T} \\ &= \sigma^{2} (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T} \mathbf{X} (\mathbf{X}^{T}\mathbf{X})^{-1} \\ &= \sigma^{2} (\mathbf{X}^{T}\mathbf{X})^{-1} \end{align*}\]

또한 \(\hat{\boldsymbol{\beta}}\) 은 반응변수 벡터 \(\mathbf{Y}\) 의 선형결합으로 표시되기 때문에 \(\mathbf{Y}\) 와 같은 정규분포를 하게 된다. 따라서 추정량 벡터의 표본분포는 다음과 같이 주어진다.

\[\begin{equation} \hat{\boldsymbol{\beta}} \sim N\left(\boldsymbol{\beta},~\sigma^{2} (\mathbf{X}^{T}\mathbf{X})^{-1}\right) \tag{4.15} \end{equation}\]

또한 개별 추정량의 분포는 다음과 같이 주어진다.

\[\begin{equation} \hat{\beta}_{j} \sim N(\beta_{j},~c_{j+1,j+1}~\sigma^{2}),~~j=0,1,\ldots,k \tag{4.16} \end{equation}\] 단, \(c_{i, j}\)은 \((\mathbf{X}^{T}\mathbf{X})^{-1}\) 행렬의 \(i\) 번째 행, \(j\) 번째 열 원소를 나타낸다.

식 (4.15)에 주어진 추정량 벡터의 표본분포를 이용해서 추론을 실시하기 위해서는 오차항 분산인 \(\sigma^{2}\) 에 대한 추정값이 반드시 필요하다. 2.3절에서 단순회귀모형의 경우에 적용되는 오차항 분산의 추정량 \(MSE\) 가 식 (2.17)으로 유도되는 것을 살펴보았다. 단순회귀모형의 경우 잔차제곱합 \(RSS\)의 자유도가 \((n-2)\) 가 된 이유는 \(\widehat{Y}_{i}\) 을 얻기 위해 추정해야할 모수의 개수가 2개이기 때문이다. 다중회귀모형에서는 추정해야 할 모수가 \(\beta_{0}, \beta_{1}, \ldots, \beta_{k}\) 의 \((k+1)\) 개가 있고, 따라서 \(RSS\)의 자유도는 \((n-k-1)\) 이 되어서 \(\sigma^{2}\) 의 불편추정량 \(MSE\) 는 다음과 같이 주어진다.

\[\begin{equation} \hat{\sigma}^{2} = \frac{RSS}{n-k-1} = MSE \tag{4.17} \end{equation}\]

식 (4.17)의 \(MSE\) 는 \(k=1\) 인 단순회귀모형의 경우도 포괄하는 일반적인 정의가 된다.

4.2.2 회귀모형의 유의성 검정

식 (4.1)의 회귀모형을 반응변수의 변동을 설명할 수 있을 것으로 생각한 모든 설명변수로 구성한 첫 번째 모형이라고 해 보자. 그러면 우리가 가장 먼저 실행할 수 있는 검정은 식 (4.1)의 회귀모형에 포함된 설명변수 중 반응변수의 변동을 설명하는데 유의적인 변수가 적어도 하나라도 있는지 알아보는 것이 될 것이며, 이것을 회귀모형의 유의성 검정이라고 한다.

귀무가설은 다음과 같고,

\[\begin{equation} H_{0}:\beta_{1}=\beta_{2}=\cdots=\beta_{k}=0 \tag{4.18} \end{equation}\]

대립가설은 다음과 같다.

\[\begin{equation} H_{1}: \text{at least one of }\beta_{j} \ne 0 \end{equation}\]

가설에 대한 검정통계량은 다음과 같다.

\[\begin{equation} F = \frac{(TSS-RSS)/k}{RSS/(n-k-1)} \end{equation}\]

귀무가설이 사실일 때 검정통게량의 분포는 \(F_{k,~n-k-1}\) 이다.

검정통계량의 구성을 살펴보자. \(TSS=\sum(y_{i}-\overline{y})^{2}\) 는 반응변수의 총변량으로서 설명변수와는 관련이 없는 변량이다. \(RSS = \sum(y_{i}-\hat{y}_{i})^{2}\) 는 잔차제곱합으로서 \(y_{i}\) 를 \(\hat{y}_{i}\) 으로 예측함으로써 발생한 오차이며, 회귀모형으로 설명이 안 된 변량이라고 할 수 있다. 따라서 총변량에서 설명이 안 된 변량을 제외한 \((TSS-RSS)\) 는 회귀모형으로 설명된 변량이라고 할 수 있다. 따라서 검정통계량은 회귀모형에 의해서 설명이 가능한 변량과 설명이 안 되는 변량의 상대적 크기라고 할 수 있으며, 검정통계량의 값이 충분히 크다면 회귀모형으로 반응변수의 변량을 설명할 수 있다는 것을 의미하기 때문에 귀무가설을 기각할 수 있게 된다.

검정통계량의 분모와 분자의 특성을 조금 더 자세하게 살펴보자. 검정통계량의 분모는 \(RSS\)를 자신의 자유도 \((n-k-1)\) 로 나눈 것으로서, 이 통계량은 회귀모형의 가정이 만족이 되면 귀무가설의 사실 여부와 관계 없이 오차항의 분산인 \(\sigma^{2}\)의 불편추정량이 된다. 즉, 가설의 사실 여부와 관계 없이 다음은 성립된다.

\[\begin{equation} E(RSS/(n-k-1))=\sigma^{2} \tag{4.19} \end{equation}\]

이것은 \(RSS/\sigma^{2}\) 이 자유도가 \((n-k-1)\) 인 \(\chi^{2}\) 분포를 하기 때문에 성립된다.

반면에, 검정통계량의 분자는 \((TSS-RSS)\) 을 자신의 자유도인 \(k\) 로 나눈 것으로서, 이 통계량은 \(H_{0}\) 이 사실인 경우에는 \(\sigma^{2}\) 의 불편추정량이 되는데, 이것은 \(H_{0}\)이 사실인 경우에 \((TSS-RSS)/\sigma^{2}\) 이 자유도가 \(k\) 인 \(\chi^{2}\) 분포를 하기 때문이다. 그러나 \(H_{0}\) 이 사실이 아닌 경우에는 \(\sigma^{2}\) 을 과대 추정하는 특성이 있다. 따라서 \(H_{0}\)이 사실이 아니면 다음과 같이 된다.

\[\begin{equation} E((TSS-RSS)/k) > \sigma^{2} \tag{4.20} \end{equation}\]

식 (4.19)와 식 (4.20)에 의해서 만일 검정통계량이 1보다 상당히 큰 값이 되면 \(H_{0}\)이 사실이 아닐 가능성이 높게 되는 것을 알 수 있다. \(H_{0}\)이 사실인 경우 검정통계량은 \(\chi^{2}\) 분포를 하는 두 통계량의 비율이 되기 때문에 \(F\) 분포를 하게 되며, 자유도는 \((k, ~n-k-1)\) 를 갖게 된다.

가설 (4.18)이 기각되면 모형에 포함된 설명변수 중 유의한 변수가 있다는 것이므로, 개별 회귀계수에 대한 유의성 검정을 진행할 수 있다. 그러나 만일 가설을 기각할 수 없다면 식 (4.1)의 회귀모형으로는 반응변수의 변동을 제대로 설명할 수 없다는 것이 된다. 이런 결과가 나오면, 다른 설명변수의 조합을 찾아 보거나, 또는 다른 형태의 회귀모형을 시도해 봐야 한다.

\(\bullet\) 자유도

2.3절에서 살펴본 제곱합의 자유도를 구하는 두 가지 방법을 적용해서 \(RSS\) 와 \((TSS-RSS)\) 의 자유도를 구해보자. 먼저 \(RSS=\sum(y_{i}-\widehat{y}_{i})^{2}\) 의 경우에는 \(n\) 개의 \(y_{i}\) 가 수집되지만 \(\widehat{y}_{i}\) 를 얻기 위해 \((k+1)\) 개 모수 \(\beta_{0}, \beta_{1}, \ldots, \beta_{k}\) 에 대한 추정이 이루어져 \((k+1)\) 개의 자유도를 잃게 되어서, 자유도는 \((n-k-1)\) 이 된다.

\((TSS-RSS)\) 의 경우에는 제곱합을 수식으로 표현해야 자유도를 구할 수 있는데, 먼저 \(TSS=\sum(y_{i}-\overline{y})^{2}\) 를 다음과 같이 분할해 보자.

\[\begin{align} \sum_{i=1}^{n}(y_{i}-\overline{y}_{i})^{2} &= \sum_{i=1}^{n}(y_{i}-\widehat{y}_{i}+\widehat{y}_{i}-\overline{y})^{2} \\ &= \sum_{i=1}^{n}(\widehat{y}_{i}-\overline{y})^{2} + \sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})^{2} +2\sum_{i=1}^{n}(\widehat{y}_{i}-\overline{y})(y_{i}-\widehat{y}_{i}) \\ &= \sum_{i=1}^{n}(\widehat{y}_{i}-\overline{y})^{2} + \sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})^{2} \tag{4.21} \end{align}\]

마지막 수식은 \(\sum_{i=1}^{n}(\widehat{y}_{i}-\overline{y})(y_{i}-\widehat{y}_{i})=0\) 이기 때문에 성립된다.

식 (4.21)의 결과에서 \((TSS-RSS)\) 는 다음과 같이 표현된다.

\[\begin{equation} TSS-RSS = \sum_{i=1}^{n}(\widehat{y}_{i}-\overline{y})^{2} \end{equation}\]

이제 수식으로 표현된 제곱합의 자유도를 구해보자. \((TSS-RSS)\) 에는 \(n\) 개의 \(\widehat{y}_{i}\) 이 있지만, \(y_{i}\) 의 경우와는 다르게 \(\widehat{y}_{i}\) 은 개별적으로 관측되는 것이 아니라 \((k+1)\) 개의 회귀계수 추정값인 \(\hat{\beta}_{0}, \hat{\beta}_{1}, \ldots, \hat{\beta}_{k}\) 에 의해서 결정되는 것이 때문에 \((k+1)\) 개의 자유도로 시작한다. 또한 \(\sum(\widehat{y}_{i}-\overline{y})=0\) 이 성립하기 때문에 하나의 제약조건이 존재하게 되어서, 자유도는 \((k+1)-1=k\) 이 된다.

4.2.3 개별 회귀계수 유의성 검정

회귀모형의 유의성 검정에서 적어도 하나의 유의적인 변수가 모형에 있다는 결론이 나면, 어떤 변수가 유의적 변수인지 확인해야 할 것이다. 개별회귀계수 유의성 검정은 식 (4.1)에 정의된 회귀모형을 구성하고 있는 각 설명변수의 유의성을 검정하는 절차이다.

가설은 다음과 같다.

\[\begin{equation} H_{0}: \beta_{j} = 0, \quad H_{1}: \beta_{j} \ne 0, \quad j = 1, \ldots, k \tag{4.22} \end{equation}\]

검정통계량은 개별 회귀계수의 추정량 \(\hat{\beta}_{j}\) 를 추정량의 표준오차로 나눈 것이다.

\[\begin{equation} t = \frac{\hat{\beta}_{j}}{SE(\hat{\beta}_{j})} \end{equation}\]

단, \(SE(\hat{\beta}_{j}) = \sqrt{c_{j+1,j+1}~MSE}\) 이며, \(c_{j+1, j+1}\) 은 (4.16)에 정의되어 있다. 귀무가설이 사실일 때 검정통계량의 분포는 \(t_{n-k-1}\) 이다.

개별 회귀계수에 대한 검정은 양측검정이 되며, 따라서 회귀계수의 신뢰구간과 밀접한 관련이 있다. 개별 회귀계수에 대한 95% 신뢰구간은 회귀계수와 추정량과 추정량의 표준오차를 이용하여 대략적으로 다음과 같이 표시된다.

\[\begin{equation} \hat{\beta}_{j} \pm 2 \cdot SE(\hat{\beta}_{j}) \end{equation}\]

만일 주어진 자료를 근거로 계산된 \(\beta_{j}\) 의 95% 신뢰구간에 0이 포함되어 있다면, 같은 자료로 계산한 검정통계량의 값은 5% 유의수준으로 구성된 기각역에 들어갈 수 없게 되어서 귀무가설을 기각할 수 없게 된다.

4.2.4 두 회귀모형의 비교

반응변수의 변동을 설명하는데 식 (4.1) 회귀모형의 설명변수를 모두 사용하는 것이 좋은지 아니면 일부분만을 사용하는 것이 더 좋은지 비교하는 절차이다. 비교하는 두 모형은 확장모형( \(\Omega\) )인 모형 (4.1)과 특정 \(q\) 개의 회귀계수에 대한 다음의 가설이 사실인 축소모형( \(\omega\) )이다.

\[\begin{equation} H_{0}: \beta_{k-q+1}=\beta_{k-q+2}=\cdots=\beta_{k}=0 \tag{4.23} \end{equation}\]

두 모형의 설명력 차이가 크지 않다면 ’모수절약의 원칙’에 의하여 축소모형을 선택하는 것이 더 좋을 것이다. 모형의 설명력 비교는 잔차제곱합을 이용할 수 있는데, 만일 확장모형의 잔차제곱합, \(RSS_{\Omega}\) 와 축소모형의 잔차제곱합, \(RSS_{\omega}\) 의 차이가 크지 않다면, 축소모형의 설명력이 확장모형 만큼 좋다는 의미가 된다.

검정통계량은 다음과 같이 구성된다.

\[\begin{equation} F=\frac{(RSS_{\omega}-RSS_{\Omega})/q}{RSS_{\Omega}/(n-k-1)} \end{equation}\]

귀무가설이 사실일 때 검정통게량의 분포는 \(F_{q,~n-k-1}\) 이다.

\(\bullet\) 예제: mtcars

4.1.1절에서 다루었던 mtcars 자료에 대한 분석을 다시 실행해 보자. 반응변수는 mpg이고, 설명변수로는 gear와 carb를 제외한 나머지 변수를 사용하며, cyl, vs, am은 요인으로 변경하자.

mtcars_4 <- mtcars |> 
  select(!c(gear, carb)) |> 
  mutate(across(c(cyl, vs, am), as.factor))

mtcars_4를 사용해서 회귀모형 fit_m3를 적합해 보자.

fit_m3 <- lm(mpg ~ ., mtcars_4)

이제 lm()으로 생성된 객체 fit_m3를 대상으로 회귀모형 추론을 위한 몇 가지 함수를 적용해 보자. 먼저 가장 빈번하게 사용되는 함수인 summary()의 결과를 확인해 보자.

summary(fit_m3)
## 
## Call:
## lm(formula = mpg ~ ., data = mtcars_4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9978 -1.3551 -0.3108  1.1992  4.1102 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 19.865583  14.801183   1.342   0.1932  
## cyl6        -1.458247   1.983190  -0.735   0.4699  
## cyl8         0.484450   3.910064   0.124   0.9025  
## disp         0.006688   0.013512   0.495   0.6255  
## hp          -0.029141   0.017182  -1.696   0.1040  
## drat         0.588059   1.503111   0.391   0.6994  
## wt          -3.155246   1.420235  -2.222   0.0369 *
## qsec         0.523235   0.690130   0.758   0.4564  
## vs1          1.237800   2.106056   0.588   0.5627  
## am1          3.000910   1.853400   1.619   0.1197  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.514 on 22 degrees of freedom
## Multiple R-squared:  0.8765, Adjusted R-squared:  0.826 
## F-statistic: 17.35 on 9 and 22 DF,  p-value: 4.814e-08

함수 lm()으로 생성된 객체에 함수 summary()를 적용시켜 얻은 결과물은 SAS나 SPSS에서 볼 수 있는 결과물과는 형식에서 차이가 있으나 많은 정보를 매우 효과적으로 보여주는 방식이라고 할 수 있다. 결과물을 하나씩 살펴보자.

먼저 Residuals:에는 잔차의 분포를 엿볼 수 있는 요약통계가 계산되어 있다. 가정이 만족된다면 잔차는 평균이 0인 정규분포를 보이게 되는데, 잔차의 요약 통계 결과 값으로 대략적인 판단을 할 수 있다.

Coefficients:에는 회귀계수의 추정값과 표준오차가 계산되어 있다. 또한 개별 회귀계수의 유의성 검정인 \(H_{0}:\beta_{j} = 0, \quad H_{1}: \beta_{j} \ne 0\) 에 대한 검정통계량의 값과 p-값이 계산되어 있다

Residual standard error:는 \(\sqrt{MSE}\) 이며, Multiple R-squared:는 결정계수 \(R^{2}\) 이고, Adjusted R-squared:는 수정결정계수의 값으로서, 모두 회귀모형의 평가 측도로 사용된다. 회귀모형의 평가측도에 대해서는 4.3절에서 살펴보겠다.

F-statistic:은 모든 회귀계수가 0이라는 가설, 즉 \(H_{0}:\beta_{1}=\beta_{2}=\cdots=\beta_{k}=0\) 에 대한 검정통계량의 값과 자유도, 그리고 p-값이 계산되어 있다.

함수 summary()로 얻어지는 결과물로 회귀모형에 대한 중요한 추론이 가능함을 알 수 있다. 한 가지 SAS나 SPSS에 익숙한 사용자들에게 아쉬울 수 있는 점은 아마도 회귀모형의 분산분석표가 출력되지 않는다는 것일 텐데, 이것은 함수 anova()를 사용함으로써 해결할 수 있다.

anova(fit_m3)
## Analysis of Variance Table
## 
## Response: mpg
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## cyl        2 824.78  412.39 65.2599 5.624e-10 ***
## disp       1  57.64   57.64  9.1218  0.006292 ** 
## hp         1  18.50   18.50  2.9279  0.101125    
## drat       1  11.91   11.91  1.8854  0.183553    
## wt         1  55.79   55.79  8.8281  0.007049 ** 
## qsec       1   1.52    1.52  0.2413  0.628165    
## vs         1   0.30    0.30  0.0478  0.828944    
## am         1  16.57   16.57  2.6216  0.119666    
## Residuals 22 139.02    6.32                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

함수 anova()는 두 회귀모형을 비교할 때에도 사용되는 함수이다. 비교되는 두 모형은 확장모형과 확장모형의 부분집합인 축소모형이어야 한다. 확장모형의 부분집합이라는 것은 확장모형의 일부 회귀계수가 0이 되는 모형을 의미하는 것으로써, 확장모형에 없는 설명변수가 축소모형에 포함되면 정상적인 비교가 이루어질 수 없다.

mtcars_4 자료에서 모든 설명변수가 포함된 확장모형 fit_m3와 변수 cyl, vs, drat를 제외한 축소모형 fit_m4를 비교해 보자.

fit_m4 <- update(fit_m3, . ~ . - cyl - vs - drat)

축소모형 fit_m4는 확장모형 fit_m3의 적합 내용을 변경한 것인데, 이런 경우에 유용하게 사용되는 함수가 update()이며, 사용법은 다음과 같다.

update(object, formula)

object는 lm()로 생성된 객체이고, formula는 object의 적합 내용 중 변경되는 모형의 공식을 지정하는 것이다. 모형 fit_m4에 적용된 모형공식은 . ~ . - cyl - vs - drat인데, 이것은 fit_m3에 적용된 모형공식 중 설명변수에서 cyl, vs, drat를 제거하는 것을 의미한다.

두 회귀모형 fit_m3와 fit_m4의 비교를 함수 anova()로 진행하자. 함수 anova()에 의한 비교 방법은 anova(축소모형, 확장모형)으로 지정되어야 한다.

anova(fit_m4, fit_m3)
## Analysis of Variance Table
## 
## Model 1: mpg ~ disp + hp + wt + qsec + am
## Model 2: mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     26 153.44                           
## 2     22 139.02  4    14.415 0.5703 0.6869

함수 anova()로 검정이 실시된 귀무가설은 변수 cyl, vs, drat의 회귀계수가 0이라는 것인데, p-값이 0.687으로 계산되어 귀무가설을 기각할 수 없게 되었다. 두 모형의 설명력에는 차이가 없다는 것이므로, 모수 절약의 원칙에 따라 축소모형을 선택하는 것이 더 좋다고 하겠다.

다중회귀모형에서 회귀계수의 유의성은 모형에 포함된 설명변수의 조합에 따라 달라질 수 있음에 유의해야 한다. 모형 fit_m3가 아닌 fit_m4를 선택했다는 것이 제외된 변수인 cyl, vs, drat가 항상 비유의적인 변수라는 것을 의미하는 것은 아니다.

lm(mpg ~ drat, mtcars_4) |> 
  summary()
## 
## Call:
## lm(formula = mpg ~ drat, data = mtcars_4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.0775 -2.6803 -0.2095  2.2976  9.0225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.525      5.477  -1.374     0.18    
## drat           7.678      1.507   5.096 1.78e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.485 on 30 degrees of freedom
## Multiple R-squared:  0.464,  Adjusted R-squared:  0.4461 
## F-statistic: 25.97 on 1 and 30 DF,  p-value: 1.776e-05

변수 drat만 사용한 회귀모형에서 drat는 유의적인 변수이고,

lm(mpg ~ drat + cyl, mtcars_4) |> 
  summary()
## 
## Call:
## lm(formula = mpg ~ drat + cyl, data = mtcars_4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3337 -1.8212 -0.1335  1.7382  6.9691 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.366      6.219   3.114  0.00423 ** 
## drat           1.793      1.509   1.188  0.24488    
## cyl6          -6.051      1.712  -3.535  0.00144 ** 
## cyl8         -10.055      1.810  -5.555 6.09e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.201 on 28 degrees of freedom
## Multiple R-squared:  0.7453, Adjusted R-squared:  0.718 
## F-statistic: 27.31 on 3 and 28 DF,  p-value: 1.83e-08

변수 drat와 cyl을 함께 사용한 모형에서는 drat가 비유의적인 변수가 되었음을 알 수 있다.

회귀계수의 신뢰구간은 함수 confint()로 계산할 수 있다. 95% 신뢰구간이 디폴트로 계산되며, 만일 신뢰수준을 변경하고자 한다면 옵션 level에 원하는 신뢰수준을 입력하면 된다. 객체 fit1에 포함되어 있는 회귀계수의 95% 신뢰구간은 다음과 같이 계산된다.

confint(fit_m4)
##                   2.5 %       97.5 %
## (Intercept) -5.66058661 34.384394537
## disp        -0.01055781  0.033033109
## hp          -0.05098537  0.008644273
## wt          -6.53883919 -1.629824922
## qsec         0.02963058  1.984163085
## am1          0.41638869  6.524518102

회귀계수의 95% 신뢰구간에 0이 포함됐다는 것은 5% 유의수준에서 가설 \(H_{0}:\beta_{j}=0, \quad H_{1}:\beta_{j} \ne 0\) 의 귀무가설을 기각할 수 없음을 의미한다. 따라서 95% 신뢰구간에 0이 포함된 변수 disp, hp는 5% 유의수준에서 비유적인 변수인 것이다.

\(\bullet\) 가변수 회귀모형 : 가변수와 연속형 변수의 상호작용 항 포함 여부

설명변수 중 범주형 변수가 있다면, 그 변수에 대한 가변수를 사용해야 하는데, 이 경우에 모형에 포함된 숫자형 설명변수와 반응변수의 관계가 범주형 변수의 모든 범주에서 동일하다고 가정하는 식 (4.13)의 모형을 사용할 수 있다. mtcars_4 자료에서 반응변수 mpg와 설명변수 qsec, am의 회귀모형을 다음과 같이 적합한다면, am의 두 범주에서 mpg와 qsec의 기울기가 동일하다고 가정한 것이 된다.

fit_m5 <- lm(mpg ~ qsec + am, mtcars_4)

모형 fit_m5의 적합 결과에 대한 그래프는 그림 4.5에서 볼 수 있다.

그러나 만일 숫자형 설명변수와 반응변수의 관계가 각 범주에서 다르다고 판단이 된다면, 식 (4.14)의 모형을 선택해서, 가변수와 연속형 변수의 상호작용 항을 포함시켜야 한다.

fit_m6 <- lm(mpg ~ qsec * am, mtcars_4)

모형 fit_m6의 모형 공식에 사용된 * 기호는 두 변수의 주 효과와 상호작용 효과를 모두 포함시키는 작용을 한다. 모형 fit_m6의 적합 결과에 대한 그래프는 그림 4.6에서 볼 수 있다.

두 모형 중 어떤 모형을 선택할 것인지는 모형 fit_m6에 포함된 상호작용 항의 유의성 여부로 결정할 수 있으나, 모형 수립 목적에 부합한 평가 기준을 사용하는 것이 더 좋은 방법이라고 하겠다.

summary(fit_m6)
## 
## Call:
## lm(formula = mpg ~ qsec * am, data = mtcars_4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4551 -1.4331  0.1918  2.2493  7.2773 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -9.0099     8.2179  -1.096  0.28226   
## qsec          1.4385     0.4500   3.197  0.00343 **
## am1         -14.5107    12.4812  -1.163  0.25481   
## qsec:am1      1.3214     0.7017   1.883  0.07012 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.343 on 28 degrees of freedom
## Multiple R-squared:  0.722,  Adjusted R-squared:  0.6923 
## F-statistic: 24.24 on 3 and 28 DF,  p-value: 6.129e-08

\(\bullet\) 예제 : modeldata::crickets

패키지 modeldata의 데이터 프레임 crickets에는 기온과 귀뚜라미 울음소리 횟수의 관계를 탐색하기 위해 관측한 자료가 입력되어 있다. 변수 temp는 기온, rate는 귀뚜라미의 분당 울음소리 횟수이며, species는 귀뚜라미 종류이다.

data(crickets, package = "modeldata")
str(crickets)
## tibble [31 × 3] (S3: tbl_df/tbl/data.frame)
##  $ species: Factor w/ 2 levels "O. exclamationis",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ temp   : num [1:31] 20.8 20.8 24 24 24 24 26.2 26.2 26.2 26.2 ...
##  $ rate   : num [1:31] 67.9 65.1 77.3 78.7 79.4 80.4 85.8 86.6 87.5 89.1 ...

세 변수의 관계 탐색을 위한 그래프를 작성해 보자. 변수 temp와 rate의 단순 산점도와 species의 효과를 살펴볼 수 있는 산점도를 모두 작성해 보자.

library(patchwork)
p1 <- ggplot(crickets, aes(x = temp, y = rate)) +
  geom_point()
p2 <- ggplot(crickets, aes(x = temp, y = rate, color = species)) +
  geom_point()
p1 + p2

그림 4.8: crickets 자료의 세 변수 관계 탐색

변수 temp와 rate 사이에는 명확한 양의 관계가 있는데, species의 두 범주 사이에는 효과 차이가 있는 것으로 보인다. 가변수를 사용한 선형회귀모형을 사용하는 것이 적절한 것으로 보이는데, 두 범주에서 변수 temp와 rate의 관계가 동일한지 여부를 확인하기 위해서 그림 4.8의 오른쪽 그래프에 회귀직선을 추가해 보자.

ggplot(crickets, aes(x = temp, y = rate, color = species)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Temperature (C)", y = "Chirp Rate (per minute)")

그림 4.9: crickets 자료의 세 변수 관계 탐색

그림 4.9에서 두 직선의 기울기는 큰 차이가 없는 것으로 보인다. 따라서 temp와 species의 상호작용 효과는 필요 없는 것으로 보이는데, 이것을 검정을 통해 확인해 보자.

식 (4.13)의 모형으로 fit.main을 적합하고, 식 (4.14)의 모형으로 fit.inter을 적합해 보자.

fit.main <- lm(rate ~ temp + species, data = crickets)
fit.inter <- lm(rate ~ temp * species, data = crickets)

모형 fit.main은 fit.inter에서 상호작용 효과만 제외된 모형으므로 fit.inter의 축소모형이 된다. 따라서 상호작용 효과의 유의성 여부는 함수 anova()를 사용해서 두 모형을 비교해서 획인할 수 있다.

anova(fit.main, fit.inter)
## Analysis of Variance Table
## 
## Model 1: rate ~ temp + species
## Model 2: rate ~ temp * species
##   Res.Df    RSS Df Sum of Sq     F Pr(>F)
## 1     28 89.350                          
## 2     27 85.074  1    4.2758 1.357 0.2542

검정 결과 p-값이 0.2542로 계산되어서, 두 모형의 설명력에는 차이가 없는 것으로 볼 수 있다. 따라서 축소모형인 fit.main을 선택하게 되고, 상호작용 효과는 비유의적이라고 결론내릴 수 있다.

물론 모형 fit.inter의 개별회귀계수 검정으로도 동일한 가설을 검정할 수 있다.

summary(fit.inter)$coef
##                          Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)           -11.0408481  4.1514800 -2.6594969 1.300079e-02
## temp                    3.7514472  0.1601220 23.4286850 1.780831e-19
## speciesO. niveus       -4.3484072  4.9616805 -0.8763981 3.885447e-01
## temp:speciesO. niveus  -0.2339856  0.2008622 -1.1649059 2.542464e-01

모형 fit.main의 적합 결과를 확인해 보자.

summary(fit.main)
## 
## Call:
## lm(formula = rate ~ temp + species, data = crickets)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0128 -1.1296 -0.3912  0.9650  3.7800 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -7.21091    2.55094  -2.827  0.00858 ** 
## temp               3.60275    0.09729  37.032  < 2e-16 ***
## speciesO. niveus -10.06529    0.73526 -13.689 6.27e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.786 on 28 degrees of freedom
## Multiple R-squared:  0.9896, Adjusted R-squared:  0.9888 
## F-statistic:  1331 on 2 and 28 DF,  p-value: < 2.2e-16

4.3 변수선택

반응변수 \(Y\) 에 대한 다음의 회귀모형에는 모든 가능한 설명변수가 다 포함되어 있다고 가정하고, 그런 의미에서 완전모형 또는 full model이라고 하자.

\[\begin{equation} Y = \beta_{0} + \beta_{1}X_{1} + \cdots + \beta_{k}X_{k} + \varepsilon \tag{4.24} \end{equation}\]

변수선택이란 반응변수의 변동을 설명할 수 있는 많은 설명변수들 중에서 ’최적’의 조합을 찾는 절차를 의미한다. 변수선택이 필요한 이유로는 모형에 불필요하게 많은 설명변수가 포함이 되면 추정 및 예측 결과에 불확실성이 증가하는 문제가 발생하기 때문이다. 정밀한 추정 및 예측을 하기 위해서는 가능한 모형을 단순하게 만들어야 한다. 또한 많은 설명변수를 포함시키다 보면 비슷한 성질의 변수가 함께 들어갈 가능성이 있는데, 그렇게 되면 다중공선성의 문제가 발생할 수 있게 된다. 적합된 모형에 대한 해석을 쉽게 하기 위해서도 가능한 단순한 모형이 필요하다고 하겠다.

변수선택 방법은 크게 세 가지로 구분할 수 있다. 첫 번째는 검정에 의하여 단계적으로 변수를 선택하는 방법이고, 두 번째는 회귀모형의 평가 측도를 근거로 변수를 선택하는 방법이며, 세 번째는 shrinkage 방법이다.

검정에 의한 변수선택은 SAS나 SPSS 등에서 일반적으로 이루어지는 방법으로서 후진소거법, 전진선택법과 단계별 선택법이 있다. 변수 선택 과정에서 요구되는 계산이 방대하지 않기 때문에 대규모의 설명변수가 있는 경우 손쉽게 중요 변수를 선택할 수 있다는 장점이 있는 방법이지만 변수의 선택과 제거가 ‘한 번에 하나씩’ 이루어지기 때문에 이른바 ‘최적’ 모형을 놓치는 경우가 발생할 수도 있다. 또한 각 단계마다 다중검정이 이루어지기 때문에 일종 오류를 범할 확률이 증가하는 검정의 정당성 문제가 발생할 수도 있다. 그리고 모형의 수립 목적이 예측인 상황에서는 변수선택 과정이 목적과 어울리지 않는다는 문제도 지니고 있는 등 많은 문제점이 있기 때문에 최근에는 거의 사용되지 않는 방법이다. 따라서 이 책에서는 두 번째와 세 번째 변수선택 방법만을 살펴보겠다.

두 번째와 세 번째 변수선택 방법을 알아보기 전에 회귀모형의 평가 측도에 대해 먼저 살펴보도록 하자.

4.3.1 회귀모형의 평가 측도

회귀모형을 적합하는 데는 나름의 목적이 있기 마련이다. 따라서 적합된 회귀모형이 그 목적에 부합한지를 적절하게 평가할 수 있는 측도가 필요하다.

결정계수( \(R^{2}\) )

결정계수는 반응변수의 총변량 중 회귀모형으로 설명된 변량의 비율을 의미한다.

\[\begin{equation} R^{2} = \frac{TSS-RSS}{TSS} = 1-\frac{RSS}{TSS} \end{equation}\]

단순회귀모형에서 결정계수는 \(Y\) 와 \(X\) 의 상관계수의 제곱과 같고, 다중회귀모형에서는 \(Y\) 와 \(\widehat{Y}\) 의 상관계수의 제곱과 같다.

반응변수의 변량을 충분히 설명하는 것이 주된 목적이라면 결정계수를 모형의 평가 측도로 사용하는 것이 좋을 것이다. 다만, 결정계수는 회귀모형에 설명변수가 추가되면 무조건 증가하는 특성이 있기 때문에 설명변수의 개수가 서로 다른 회귀모형을 비교하는 평가 측도로는 적절하지 않으며, 반드시 설명변수의 개수가 같은 회귀모형을 비교할 때 사용해야 한다.

수정결정계수(\(adj~R^{2}\) )

회귀모형에 설명변수가 추가되면 무조건 값이 증가하는 결정계수의 문제를 보완한 측도이다. 추가된 설명변수가 모형의 설명력 향상에 도움이 되는 경우에만 값이 증가되도록 결정계수를 수정한 측도이다.

\[\begin{equation} adj~R^{2}=1-\frac{RSS/(n-k-1)}{TSS/(n-1)} = 1-\left( \frac{n-1}{n-k-1} \right)(1-R^{2}) \end{equation}\]

Residual standard error(RSE)

잔차제곱합, \(RSS\) 는 반응변수의 총변량 중에 회귀모형으로 설명이 되지 않는 변량으로서, 결정계수의 값을 결정하는 요소이며, 회귀모형에 설명변수가 추가되면 무조건 감소하는 특성을 가지고 있다. 이러한 특성은 다음과 같이 잔차제곱합을 자신의 자유도로 나누면 수정이 된다.

\[\begin{equation} RSE = \sqrt{\frac{RSS}{n-k-1}} \tag{4.25} \end{equation}\]

수정결정계수와 실질적으로 동일한 특성을 가지고 있는 측도로서, 추가된 설명변수가 모형의 설명력 향상에 도움이 되는 경우에만 값이 감소하는 특성을 가지고 있다.

Information criteria

회귀모형에서 오차항이 정규분포를 한다는 가정을 하고 있기 때문에 Likelihood function을 구성할 수 있는데, 회귀모형에 설명변수를 추가함으로써 모수의 개수가 증가하게 되면 Maximum likelihood의 값도 자연스럽게 증가하게 된다. 이것은 결정계수의 특성과 유사한 것이며, 따라서 설명변수 증가에 대한 penalty를 부과하는 측도를 생각해 볼 수 있다.

AIC와 BIC는 다음과 같이 정의된다.

\[\begin{align*} AIC &= -2\log \hat{L} + 2K \\ BIC &= -2\log \hat{L} + k\log n \end{align*}\]

단, \(\hat{L}\) 은 maximum likelihood function이고, \(n\) 은 자료의 크기, \(k\) 는 설명변수의 개수이다. 회귀모형에 설명변수가 추가되면, \(\hat{L}\) 의 값은 증가하게 되고, 따라서 \(-2\log \hat{L}\) 은 감소한다. 설명변수 증가에 대한 penalty로써 AIC는 \(2k\)를, BIC는 \(k \log n\)를 사용하고 있는데, 설명변수를 추가함으로써 증가된 penalty의 값보다 모형의 적합도 향상으로 감소된 \(-2 \log \hat{L}\) 의 값이 더 크다면 AIC나 BIC는 감소하게 된다. 따라서 추가된 설명변수가 모형의 적합도 향상에 도움이 되는 경우에만 값이 감소하는 특성을 가지고 있다.

\(C_{p}\) 통계량

정확한 예측 결과가 중요한 평가 기준이라면 다음에 주어진 예측오차제곱합을 기준으로 사용할 수 있을 것이며,

\[\begin{equation} \frac{1}{\sigma^{2}}\sum_{i} E\left(\widehat{Y}_{i}-E(Y_{i})\right)^{2} \end{equation}\]

추정은 \(C_{p}\) 통계량으로 할 수 있다.

\[\begin{equation} C_{p} = \frac{RSS_{p}}{\hat{\sigma}^{2}} + 2(p+1) - n \end{equation}\]

단, \(\hat{\sigma}^{2}\) 는 식 (4.24)의 완전모형에서 추정한 오차항의 분산이며, \(RSS_{p}\) 는 설명변수의 개수가 \(p\) 개인 모형의 잔차제곱합이다. 설명변수의 개수가 \(p\) 개인 모형의 경우에는 \(E(RSS_{p})=(n-p-1)\sigma^{2}\)가 되기 때문에, \(E(C_{p}) \approx (p+1)\) 가 될 것이다. 만일 주어진 회귀모형이 자료를 잘 설명하지 못하는 모형이라면, 잔차제곱합이 매우 큰 값이 될 것이 때문에 \(C_{p}\) 값이 \((p+1)\) 보다 큰 값이 될 것이다. 따라서 가능한 작은 \(p\) 값에서 \(C_{p}\) 의 값이 \((p+1)\) 과 비슷하거나 작게 되는 모형을 선택하면 된다. 완전모형에서는 \(C_{p}=p+1\) 이 된다.

\(\bullet\) Cross-validation(CV)

CV는 회귀모형의 예측력을 평가할 수 있는 재표본추출(resampling) 기법이다. 지금까지 살펴본 모형 평가 측도는 적합된 모형을 대상으로 계산하는 것인데, 문제는 어떤 자료를 사용해서 모형을 적합하고, 평가 측도값을 계산하는 것이 가장 적절한 것인지이다.

예측 목적으로 설정된 회귀모형의 경우에는 모형 적합에 사용된 자료의 변동 설명보다 적합에 사용되지 않은 새로운 자료에 대한 예측 오차가 더 중요한 평가 요소라고 할 수 있다. 통계모형 적합에 사용되지 않은 새로운 자료에 대한 예측 오차를 test error라고 하는데, 회귀모형의 예측 결과에 대한 정당성을 확보하기 위해서는 낮은 test error가 필수적이다.

CV는 모형 적합에 사용되는 자료인 training data를 이용해서 test error를 추정하기 위한 재표본추출 기법이다. 몇 가지 조금 다른 방법이 있는데, 여기에서는 k-fold CV라는 방법에 대해서만 살펴보겠다. k-fold CV는 training data를 비슷한 크기의 k개 그룹(fold)으로 분리하는 것으로 시작한다. 분리된 k개의 그룹 중 첫 번째 그룹의 자료를 제외하고, 나머지 (k-1)개 그룹의 자료를 사용해서 적합한 회귀모형으로 적합 과정에서 제외된 첫 번째 그룹에 대한 예측을 실시하고, 적절한 평가 측도를 사용해서 예측 오차를 계산한다. 이어서 각 그룹의 자료를 하나씩 차례로 제외하고 나머지 자료로 적합 및 예측하는 과정을 반복하여 k번의 예측을 실시한다. Test error는 k번의 예측에서 각각 발생된 예측 오차의 평균으로 추정하게 된다.

4.3.2 평가 측도에 의한 변수선택

평가 측도에 의한 방법은 4.3.1절에서 살펴본 \(AIC\), \(BIC\), \(adj~R^{2}\) 및 \(C_{p}\) 통계량 등을 기반으로 ‘최적’ 변수를 선택하는 방법이다. 모형에 포함시킬 수 있는 설명변수가 \(k\) 개 있다면, 적합 가능한 모형의 개수는 \(2^{k}\) 개가 된다. 이 경우에 \(2^{k}\) 개의 모형을 다 적합시키고 특정 평가 측도를 기준으로 최적 모형을 선택하는 이른바 ‘Best subset selection’ 방법을 적용해 볼 수 있다. 하지만, 이 방법은 \(k\) 의 값이 커진다면 현실적으로 적용하기 쉽지 않다는 문제를 갖고 있다.

다른 방법으로는 특정 평가 측도를 기준으로 변수를 하나씩 모형에 추가하거나 제거하는 ‘Stepwise selection’ 방법을 생각해 볼 수 있다. 이 방법은 검정에 의한 변수선택에서 적용되는 방법으로서 ‘최적’ 모형을 찾지 못할 가능성이 있는 방법이다. 나름의 장점과 단점이 존재하는 방법 중 주어진 상황에 가장 적합한 방법을 선택해서 사용해야 할 것이다.

\(\bullet\) Best subset selection

\(k\) 개 설명변수의 모든 가능한 조합 중 모형의 성능을 최대로 할 수 있는 조합을 찾는 방법이다. 일반적으로 적용되는 방식은 다음과 같다.

\(i=1,2,\ldots,k\) 에 대하여 \(i\) 개 설명변수가 있는 \(\binom{~k~}{i}\) 개 모형을 적합하고, 그 중 결정계수의 값이 가장 큰 모형을 \(M_{i}\) 로 지정
\(M_{0}, M_{1}, \ldots, M_{k}\) 를 대상으로 \(AIC\), \(BIC\), \(C_{p}\), \(adj~R^{2}\) 등을 기준으로 ‘최적’ 모형 선택

\(\bullet\) Stepwise selection

‘한 번에 하나씩’ 이루어지는 선택 과정으로 인하여 ‘최적’ 변수의 조합을 찾지 못할 수 있다는 문제가 있지만, 설명변수의 개수가 많은 경우에 비교적 빠르게 변수선택을 할 수 있다는 장점이 있는 방법이다. 적용할 수 있는 방법에는 모형에 변수를 하나씩 추가하는 Forward stepwise selection과 모형에서 하나씩 변수를 제거하는 Backward elimination, 그리고 두 가지 방식이 혼합된 Hybrid stepwise selection이 있다.

1. Forward stepwise selection

절편만 있는 모형 \(M_{0}~(Y=\beta_{0}+\varepsilon)\) 에서 시작하여 단계적으로 변수를 하나씩 모형에 추가하는 방법이다. 추가할 변수를 선택하는 절차는 다음과 같이 두 가지 방식이 있다.

\(\bullet\) 방식 1

\(AIC\) 등을 기준으로 모형의 성능을 최대로 개선시킬 수 있는 변수를 하나씩 모형에 추가
더 이상 모형의 성능을 개선시킬 변수가 없으면, 절차 종료

\(\bullet\) 방식 2

\(i=0, \ldots, k-1\) 에 대하여 \(M_{i}\) 에 하나의 설명변수를 추가한 \((k-i)\) 개 모형 중 결정계수 값이 가장 큰 모형을 \(M_{i+1}\) 으로 지정
\(M_{0}, M_{1}, \ldots, M_{k}\) 를 대상으로 \(AIC\), \(BIC\), \(C_{p}\), \(adj~R^{2}\) 등을 기준으로 ‘최적’ 모형을 하나 선택

Forward stepwise selection에서는 일단 모형에 포함된 변수에 대해서는 추가적인 확인 과정이 없다. 그러나 이러한 특성은 개별 설명변수의 영향력은 모형에 추가적으로 포함되는 다른 설명변수에 의해 변할 수 있기 때문에 문제가 될 수 있다.

2. Backward stepwise elimination

모든 설명변수가 포함된 full model, \(M_{k}\) 에서 시작하여 단계적으로 변수를 하나씩 제거하는 방법이다. 제거할 변수를 선택하는 절차는 다음과 같이 두 가지 방식이 있다.

\(\bullet\) 방식 1

제거되면 모형의 성능이 최대로 개선되는 변수를 하나씩 제거
모형에 포함된 변수 중 어떤 변수라도 제거되면 모형의 성능이 더 나빠질 때 절차 종료

\(\bullet\) 방식 2

\(i=k, k-1, \ldots, 1\) 에 대하여 \(M_{i}\) 에서 하나씩 변수가 제거된 모구 \(i\) 개 모형 중 결정계수의 값이 가장 큰 모형을 \(M_{i-1}\) 으로 지정
\(M_{0}, M_{1}, \ldots, M_{k}\) 를 대상으로 \(AIC\), \(BIC\), \(C_{p}\), \(adj~R^{2}\) 등을 기준으로 ‘최적’ 모형을 하나 선택

Backward stepwise elimination에서는 일단 제거된 변수는 다시 모형에 포함될 수 없는 방식이다. 그러나 이러한 특성은 개별 설명변수의 영향력은 모형에서 다른 설명변수가 제거되면 변할 수 있기 때문에 문제가 될 수 있다.

3. Hybrid stepwise selection

Forward stepwise selection과 Backward stepwise elimination는 일단 모형에 포함되거나 혹은 모형에서 제거된 변수에 대한 추가적인 고려가 불가능한 방식이다. 이러한 특성으로 인하여 ‘최적’ 모형을 찾지 못할 가능성이 있는데, 이 문제는 두 방식의 특성을 혼합시킨 방식을 적용하면 해결할 수 있다.

Hybrid forward selection : 각 단계별로 변수를 추가한 후에 모형에 포함되어 있는 변수를 대상으로 backward elimination 기법을 적용해서 변수 제거
Hybrid backward elimination : 각 단계별로 변수를 제거한 후에 이미 제거된 변수를 대상으로 forward selection 기법으로 모형에 추가될 변수 선택

\(\bullet\) 변수선택을 위한 함수

Best subset selection은 함수 leaps::regsubsets()로 할 수 있다. 기본적인 사용법은 다음과 같다.

regsubsets(formula, data, nbest = 1, nvmax = 8)

formula와 data는 함수 lm()에서 사용한 방식과 동일하다. 설명변수의 개수가 동일한 모형 중 결정계수의 값이 가장 큰 nbest개 모형을 선택하는데, 이렇게 선택된 모형을 대상으로 다양한 평가 측도를 기준으로 비교할 수 있다. 또한 최대 nvmax개의 설명변수가 포함된 모형을 대상으로 비교를 하기 때문에, nvmax에는 설명변수의 개수를 지정할 필요가 있다.

Stepwise selection은 함수 MASS::stepAIC()로 할 수 있다. 기본적인 사용법은 다음과 같다.

stepAIC(object, scope, direction, trace = 1, k = 2)

object는 stepwise selection을 시작하는 회귀모형을 나타내는 lm 객체를 지정한다. Forward selection의 경우에는 절편만이 있는 모형인 lm(y ~ 1, data)가 될 것이고, Backward elimination의 경우에는 모든 설명변수가 포함된 모형인 lm(y ~ ., data)가 될 것이다. scope는 탐색 범위를 나타내는 upper와 lower를 요소가 갖고 있는 리스트를 지정한다. direction은 탐색 방식을 지정하는데, "forward"는 모형에 포함된 변수는 계속 유지하면서 forward selection 방법을 진행하며, "backward"는 일단 모형에서 제거된 변수는 최종 모형에 포함시키지 않는 방식으로 backward elimination을 진행한다. 또한 "both"는 hybrid stepwise selection 방법을 진행한다. direction의 디폴트는 scope를 지정하면 "both"가 되지만, scope가 생략되면 "backward"가 된다. trace는 탐색 과정의 출력 여부를 지정하는 것으로 출력되는 것이 디폴트이며, FALSE 또는 0을 지정하면 출력되지 않는다. k는 AIC나 BIC의 계산 과정에서 모형에 포함된 변수의 개수만큼 불이익을 주기 위한 상수로서, 디폴트는 AIC 계산을 위한 k = 2이다. 만일 BIC에 의한 단계별 선택법을 적용하고자 한다면 k = log(n)으로 수정해야 한다. 단, n은 데이터 프레임의 행 개수이다.

\(\bullet\) 예제: mtcars

4.1.1절에서 다루었던 mtcars를 대상으로 변수선택을 진행해 보자. mtcars를 대상으로 앞서 이루어진 주된 분석은 다음과 같다.

mtcars_4 <- mtcars |> 
  select(!c(gear, carb)) |> 
  mutate(across(c(cyl, vs, am), as.factor))

Best subset selection 방법을 함수 leaps::regsubsets()로 진행해 보자.

library(leaps)
fits <- regsubsets(mpg ~ ., mtcars_4, nvmax = 9)

함수 regsubsets()에 요인을 설명변수로 입력하면, 각 요인에 대한 가변수를 자동으로 생성한다. 그런데 여기에서 주의할 점은 각 가변수를 개별적인 설명변수로 취급한다는 것이다.

모형 fits에 입력된 설명변수는 연속형 변수 5개와 요인 3개인데, 요인의 경우에는 각 가변수를 개별 변수로 인식하기 때문에 cyl에 대한 가변수 2개, vs와 am에 대한 가변수 각각 1개씩 모두 4개의 변수가 있는 것으로 인식한다. 따라서 모형 fits에는 9개의 설명변수가 입력되기 때문에 nvmax = 9를 지정한 것이다.

설명변수의 개수가 \(i\) 개 (\(i=1, \ldots, 9\)) 인 모형 중 결정계수가 가장 높은 nbest개의 모형은 다음과 같이 함수 summary()로 출력할 수 있다. 모형에 포함되는 변수에는 "*"가 표시되어 있다.

summary(fits)
## Subset selection object
## Call: regsubsets.formula(mpg ~ ., mtcars_4, nvmax = 9)
## 9 Variables  (and intercept)
##      Forced in Forced out
## cyl6     FALSE      FALSE
## cyl8     FALSE      FALSE
## disp     FALSE      FALSE
## hp       FALSE      FALSE
## drat     FALSE      FALSE
## wt       FALSE      FALSE
## qsec     FALSE      FALSE
## vs1      FALSE      FALSE
## am1      FALSE      FALSE
## 1 subsets of each size up to 9
## Selection Algorithm: exhaustive
##          cyl6 cyl8 disp hp  drat wt  qsec vs1 am1
## 1  ( 1 ) " "  " "  " "  " " " "  "*" " "  " " " "
## 2  ( 1 ) " "  " "  " "  "*" " "  "*" " "  " " " "
## 3  ( 1 ) " "  " "  " "  " " " "  "*" "*"  " " "*"
## 4  ( 1 ) "*"  " "  " "  "*" " "  "*" " "  " " "*"
## 5  ( 1 ) "*"  " "  " "  "*" " "  "*" " "  "*" "*"
## 6  ( 1 ) "*"  " "  " "  "*" " "  "*" "*"  "*" "*"
## 7  ( 1 ) "*"  " "  "*"  "*" " "  "*" "*"  "*" "*"
## 8  ( 1 ) "*"  " "  "*"  "*" "*"  "*" "*"  "*" "*"
## 9  ( 1 ) "*"  "*"  "*"  "*" "*"  "*" "*"  "*" "*"

설명변수가 하나인 경우에는 wt, 두 개인 경우에는 wt와 hp, 세 개인 경우에는 wt, qsec, am이 선택된 것을 알 수 있다. 2개 이상의 가변수가 사용된 범주형 변수의 경우에는 하나의 가변수만 포함되어도 해당 범주형 변수를 모형에 포함시키는 것이 적절할 것이다.

이제 선택된 9개 모형을 대상으로 평가 측도를 비교해 보자. 평가 측도의 비교는 객체 fits를 함수 plot()에 입력하면 된다. 평가 측도는 scale에 지정할 수 있는데, 가능한 키워드는 "bic", "Cp", "adjr2", "r2"이며, 디폴트는 "bic"이다.

plot(fits)

그림 4.10: BIC에 의한 비교 결과

그림 4.10의 X축에는 설명변수가 표시되어 있고, Y축에는 비교 대상이 되는 각 모형의 \(BIC\) 의 값이 표시되어 있다. 즉, 그래프의 각 행은 비교 대상이 되는 모형을 나타내는 것이며, X축에 표시된 설명변수가 해당 모형에 포함이 된 것이면 직사각형에 색이 채워진다. 위에서 첫 번째 모형이 \(BIC\) 의 값이 가장 작은 모형이며, 변수 wt, qsec, am이 포함된 모형이다.

plot(fits, scale = "adjr2")

그림 4.11: 수정결정계수에 의한 비교 결과

수정결정계수를 평가 측도로 변경하여 변수선택을 진행해 보자. 수정결정계수가 가장 큰 모형으로는 변수 hp, wt, vs, am과 cyl이 포함된 모형이 선택되었다.

plot(fits, scale = "Cp")

그림 4.12: Cp에 의한 비교 결과

\(C_{p}\) 가 평가 측도인 경우, 선택된 모형은 수정결정계수의 경우와 큰 차이가 없음을 알 수 있다.

Stepwise selection 방법을 함수 MASS::stepAIC()로 진행해 보자. 먼저 모든 설명변수가 포함된 모형과 절편만 있는 모형에 대한 lm 객체를 생성시키자.

fit_full <- lm(mpg ~ ., mtcars_4)
fit_null <- lm(mpg ~ 1, mtcars_4)

AIC에 의한 forward selection으로 변수선택을 진행해 보자. fit_null로 시작하고 scope의 upper가 fit_full로 지정하였기 때문에 hybrid forward selection이 수행된다. 변수 wt, cyl, hp, am이 선택되었다.

MASS::stepAIC(fit_null,
        scope = list(lower = fit_null, upper = fit_full),
        trace = FALSE)
## 
## Call:
## lm(formula = mpg ~ wt + cyl + hp + am, data = mtcars_4)
## 
## Coefficients:
## (Intercept)           wt         cyl6         cyl8           hp          am1  
##    33.70832     -2.49683     -3.03134     -2.16368     -0.03211      1.80921

AIC에 의한 backward elimination으로 변수선택을 진행해 보자. Backward elimination은 full model을 시작 모형으로 지정하면 실행된다. AIC에 의한 forward selection과 같은 결과를 보이고 있다.

MASS::stepAIC(fit_full, direction = "both", trace = FALSE) 
## 
## Call:
## lm(formula = mpg ~ cyl + hp + wt + am, data = mtcars_4)
## 
## Coefficients:
## (Intercept)         cyl6         cyl8           hp           wt          am1  
##    33.70832     -3.03134     -2.16368     -0.03211     -2.49683      1.80921

BIC에 의한 forward selection으로 변수선택을 진행해 보자. 변수 wt와 hp만 선택되어 AIC에 의한 방법과는 다른 결과가 나왔다.

MASS::stepAIC(fit_null, 
          scope = list(lower = fit_null, upper = fit_full),
          k = log(nrow(mtcars_4)), 
          trace = FALSE)
## 
## Call:
## lm(formula = mpg ~ wt + hp, data = mtcars_4)
## 
## Coefficients:
## (Intercept)           wt           hp  
##    37.22727     -3.87783     -0.03177

BIC에 의한 backward elimination으로 변수선택을 진행해 보자. 변수 wt, qsec, am이 선택되었다.

MASS::stepAIC(fit_full, direction = "both", 
          k = log(nrow(mtcars_4)), trace = FALSE)
## 
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars_4)
## 
## Coefficients:
## (Intercept)           wt         qsec          am1  
##       9.618       -3.917        1.226        2.936

4.3.3 Shrinkage 방법

4.3.2절에서 살펴본 변수선택 방식은 설명변수들의 부분집합으로 구성되는 다양한 회귀모형을 적합하고, 그 모형의 평가측도로 ‘최적’ 모형을 선택하는 것이다. 이 방식의 대안으로 최근 많이 사용되는 shrinkage 방법은 회귀계수의 크기에 제약 조건을 부과해서 중요하지 않은 변수의 회귀계수를 0으로 줄어들게 만드는 방법으로서 ridge 회귀모형과 lasso가 여기에 속한다.

4.3.3.1 Ridge 회귀모형

식 (4.24)에 주어진 회귀모형의 회귀계수를 최소제곱추정법으로 추정하는 것은 \(RSS\) 를 최소화하는 회귀계수를 구하는 것을 의미한다.

\[\begin{equation} RSS = \sum_{i=1}^{n} \left(y_{i}-\beta_{0}-\beta_{1}x_{1i}-\cdots-\beta_{k}x_{ki} \right)^{2} \end{equation}\]

Ridge 회귀모형에서 회귀계수를 추정하는 방법은 단순히 \(RSS\) 만을 최소화하는 것이 아니라 회귀계수의 크기를 제한하는 조건이 추가하는 것이다. Ridge 회귀모형의 회귀계수는 식 (4.26)을 최소화시키는 값이 된다.

\[\begin{equation} \sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\sum_{j=1}^{k}\beta_{j}x_{ji} \right)^{2} + \lambda \sum_{j=1}^{k}\beta_{j}^{2} \tag{4.26} \end{equation}\] 단, \(\lambda\) 는 조율모수(tuning parameter)이다.

식 (4.26)은 두 개의 성분으로 구성되어 있는데, 첫 번째 성분은 \(RSS\) 로서 자료에 잘 적합이 되면 작은 값을 갖게 되는 성분이지만, 두번째 성분은 회귀계수의 크기에 비례해서 작아지는 성분이다. 두 번째 성분을 수축 패널티(shrinkage penalty) 성분이라고 하는데, 그것은 대부분의 회귀계수 값이 \(0\) 에 가까운 값으로 수축되어야 작은 값을 갖기 때문이다.

두 번째 성분의 영향력은 조율모수 \(\lambda\) 의 크기에 의해서 조절되는데, 만일 \(\lambda = 0\) 이면 최소제곱추정법과 동일한 결과가 나올 것이고, 반대로 대단히 큰 값이 주어지면 자료의 적합 정도에 관계없이 모든 회귀계수가 \(0\) 에 가까운 값을 갖게 될 것이다. 조율모수 \(\lambda\) 의 최적 값을 구하는 것은 ridge 회귀모형에서 가장 중요한 작업이 되는데, 일반적으로 cross-validation으로 구하게 된다.

Ridge 회귀모형은 조율모수 \(\lambda\) 의 값이 커짐에 따라 자료에 대한 적합 정도는 점점 떨어지게 되는데, 이러한 특성으로 편기(bias)가 증가하게 된다. 하지만 같은 이유로 회귀계수의 추정값은 자료의 변동에 덜 민감하게 되고, 결과적으로 분산이 감소하게 된다. 따라서 최소제곱추정법에 의해 추정된 회귀계수의 분산이 매우 큰 값을 갖게 되는 경우에는 ridge 회귀모형이 대안으로 사용될 수 있다.

Ridge 회귀모형의 회귀계수 추정값은 영향력이 작은 변수의 경우에는 상대적으로 0에 더 근접한 값으로 수축되지만, 정확하게 0의 값이 되는 것은 아니다. 따라서 변수선택의 목적으로 사용하기에는 적절하지 않은 방법이 된다.

4.3.3.2 Lasso

Lasso(Least Absolute Shrinkage and Selection Operator)도 ridge 회귀모형처럼 \(RSS\) 와 회귀계수의 크기를 제한하는 조건이 추가된 변량을 최소화하는 방식으로 회귀계수를 추정한다. Lasso 회귀계수는 식 (4.27)을 최소화시키는 값이 된다.

\[\begin{equation} \sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\sum_{j=1}^{k}\beta_{j}x_{ji} \right)^{2} + \lambda \sum_{j=1}^{k}|\beta_{j}| \tag{4.27} \end{equation}\]

Lasso도 ridge 회귀모형의 경우처럼 조율모수 \(\lambda\) 의 값이 커지게 되면 회귀계수 추정값이 수축하게 된다. 다만, 두 모형의 차이점은 수축 패널티 성분의 형태인데, ridge 회귀모형에서는 \(\beta_{j}^{2}\) 을 사용하고, lasso에서는 \(|\beta_{j}|\) 을 사용한다. 제곱 형태의 패널티를 사용하면 정확하게 0으로 수축시킬 수 없지만, 절대값 형태의 패널티를 사용하면 0으로 수축시키는 것이 가능하게 된다. 따라서 lasso에서는 영향력이 작은 변수의 회귀계수가 정확하게 0의 값을 갖게 되어, 실질적인 변수선택을 실시하는 것이 된다.

\(\bullet\) Shrinkage 모형 적합을 위한 함수

Ridge와 lasso 모형은 패키지 glmnet의 함수를 사용해서 적합할 수 있다. 기본적인 적합 함수는 glmnet()이며 사용법은 다음과 같다.

glmet(x, y, alpha = 1)

x에는 설명변수의 행렬을 지정하고, y에는 반응변수의 벡터를 지정한다. alpha는 0과 1 사이의 값을 갖는 elasticnet 모수인데, alpha = 1이면 lasso, alpha = 0이면 ridge 패널티가 적용된다. 모형 적합은 조율모수 \(\lambda\) 에 여러 값을 차례로 입력해서 이루어진다.

최적 \(\lambda\) 값을 cross-validation으로 추정한 모형을 적합하기 위해서는 함수 cv.glmnet(x, y, alpha)을 사용하는 것이 더 간편할 수 있다.

\(\bullet\) 예제: mtcars

4.3.2절에서 다루었던 mtcars를 대상으로 shringkage 모형을 적합해 보자. mtcars를 대상으로 앞서 이루어진 주된 분석은 다음과 같다.

mtcars_4 <- mtcars |> 
  select(!c(gear, carb)) |> 
  mutate(across(c(cyl, vs, am), as.factor))

패키기 glmnet의 함수를 사용하기 위해 설명변수로 이루어진 행렬과 반응변수 벡터를 생성해보자.

X <- model.matrix(mpg ~ ., mtcars_4)[,-1]
Y <- mtcars_4$mpg

설명변수의 행렬은 각 변수를 함수 cbind()로 연결해서 만들 수 있지만, 변수의 개수가 많거나 요인이 포함되어 있는 경우에는 불편한 방법이 된다. 함수 model.matrix()는 회귀모형의 design matrix인 \({\bf X}\) 행렬을 생성하는 기능이 있는 함수이다. model.matrix(formula, data)의 형태로 lm()과 동일한 방식으로 모형을 정의하면 된다. 생성된 \({\bf X}\) 행렬의 첫 번째 열은 절편에 대한 것이기 때문에 설명변수의 행렬에 해당되지 않으므로 제외한다.

먼저 ridge 회귀모형을 적합해보자.

library(glmnet)
fit_R <- glmnet(X, Y, alpha = 0)

조율모수 \(\lambda\) 값의 변화에 따른 회귀계수 추정값의 변화 그래프는 fit_R을 함수 plot()에 입력하면 작성할 수 있다.

plot(fit_R, xvar = "lambda")

그림 4.13: lambda 값의 변화에 따른 ridge 회귀계수 추정값의 변화

그림 4.13은 X축에 표시된 \(\log \lambda\) 의 값 변화에 따른 회귀계수 추정 결과의 변화를 선으로 나타낸 그래프이다. \(\lambda\) 값이 커지면서 회귀계수가 0으로 수축되고 있는 것을 볼 수 있다. 그래프 위쪽의 눈금은 0이 아닌 회귀계수의 개수를 나타내고 있는데, ridge 회귀계수는 매우 큰 \(\lambda\) 값에 대해서도 모두 0이 되지는 않고 있음을 알 수 있다.

Lasso 모형을 적합하고, \(\lambda\) 의 값 변화에 따른 회귀계수 추정 결과의 변화를 나타내 보자.

fit_L <- glmnet(X, Y, alpha = 1)
plot(fit_L, xvar = "lambda")

그림 4.14: lambda 값의 변화에 따른 lasso 회귀계수 추정값의 변화

그림 4.14을 보면 \(\log \lambda\) 값이 증가함에 따라서 0이 아닌 회귀계수의 개수가 줄어들고 있음을 알 수 있다.

최적 \(\lambda\) 값을 cross-validation으로 추정한 모형을 함수 cv.glmnet()으로 적합해 보자. 먼저 ridge 회귀모형을 적합해보자.

cvfit_R <- cv.glmnet(X, Y, alpha = 0)
cvfit_R
## 
## Call:  cv.glmnet(x = X, y = Y, alpha = 0) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min  0.899    94   6.547 1.694       9
## 1se  6.345    73   8.192 2.590       9

10-fold cross-validation을 사용하는 것이 디폴트이며, 반응변수가 연속형인 경우에는 MSE를 최소화시키는 \(\lambda\) 값을 추정한다. min에 해당하는 Lambda는 CV MSE를 최소화시키는 \(\lambda\) 값이고, 1se에 해당하는 Lambda는 최소 CV MSE \(\pm~~ 1 \times SE\) 범위 안에서 최대 \(\lambda\) 값이 된다. 1se의 \(\lambda\) 값은 최소 CV MSE와는 큰 차이가 없으면서도 가능한 더 많은 회귀계수를 0으로 수축시킬 수 있는 조율모수가 된다.

\(\lambda\) 값의 변화에 따른 CV의 결과는 다음과 같이 확인할 수 있다.

plot(cvfit_R)

그림 4.15: lambda 값의 변화에 따른 CV MSE의 변화

왼쪽에서 첫 번째 수직 점선이 최소 CV MSE에 해당하는 \(\lambda\) 값이고, 두 번째 수직 점선이 1SE에 해당하는 \(\lambda\) 값이다.

이번에는 lasso 모형을 적합해 보자.

cvfit_L <- cv.glmnet(X, Y, alpha = 1)
cvfit_L
## 
## Call:  cv.glmnet(x = X, y = Y, alpha = 1) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min 0.1983    36   7.779 2.801       7
## 1se 1.1617    17  10.195 5.096       5

\(\lambda\) 값의 변화에 따른 CV의 결과도 확인해 보자.

plot(cvfit_L)

그림 4.16: lambda 값의 변화에 따른 CV MSE의 변화

최적 \(\lambda\) 값으로 추정된 회귀계수 값은 함수 coef()에 cv.glmnet 객체를 입력하면 확인할 수 있다. s = "lambda.min"을 추가하면 최소 CV MSE에 해당하는 \(\lambda\) 값으로 추정된 회귀계수가 출력되고, 생략하면 디폴트 값인 s = "lambda.1se"가 지정되어 1SE에 해당하는 \(\lambda\) 값으로 추정된 회귀계수가 출력된다.

먼저 ridge 회귀계수의 추정값을 출력해 보자.

coef(cvfit_R)
## 10 x 1 sparse Matrix of class "dgCMatrix"
##                       s1
## (Intercept) 19.022276894
## cyl6        -1.079993964
## cyl8        -1.072884386
## disp        -0.006325128
## hp          -0.013135723
## drat         1.107153854
## wt          -1.124668959
## qsec         0.208350553
## vs1          1.119832155
## am1          1.447003442

Lasso 회귀계수의 추정값을 출력해 보자.

coef(cvfit_L)
## 10 x 1 sparse Matrix of class "dgCMatrix"
##                        s1
## (Intercept) 32.0056920950
## cyl6         .           
## cyl8         .           
## disp        -0.0008627785
## hp          -0.0202233183
## drat         0.2171493960
## wt          -2.9820129233
## qsec         .           
## vs1          0.1448971143
## am1          .

회귀계수의 추정값이 점으로 찍힌 변수는 해당 회귀계수가 0이라는 것을 의미한다. Ridge 회귀계수와 lasso 회귀계수의 추정 결과는 최소제곱추정에 의한 추정결과와는 당연히 다르며, 회귀계수 추정에 대한 SE로 계산되지 않는다. 또한 cross-validation에 의한 결과이기 때문에 추정할 때마다 추정결과가 조금 다르게 된다.

4.4 연습문제

1. 반응변수 \(Y\) 에 대하여 설명변수 \(X_{1}, \ldots, X_{k}\) 의 회귀모형 \(\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\) 을 설정하고자 한다.

반응변수의 관찰값 벡터 \(\mathbf{Y}\) 를 정의하라.
설명변수의 관찰값 행렬이며, design matrix라고 불리는 \(\mathbf{X}\) 를 정의하라.
모수 벡터 \(\boldsymbol{\beta}\) 를 정의하라.
오차항 벡터 \(\boldsymbol{\varepsilon}\) 를 정의하라.
모수 벡터 \(\boldsymbol{\beta}\) 의 최소제곱추정량 벡터 \(\hat{\boldsymbol{\beta}}\) 은 무엇인가?

2. 다중회귀모형 \(y_{i}=\beta_{0}+\beta_{1}x_{1i}+\cdots+\beta_{k}x_{ki}+\varepsilon_{i}\), \(i=1,\ldots,n\)에 대한 추론 및 모형 평가를 실시하고자 한다.

회귀모형의 유의성 검정이란 무엇이며, 어떤 가설에 대한 검정인지 설명하라.
회귀모형의 유의성 검정에 대한 검정통계량은 \(F=\frac{TSS-RSS}{k}/\frac{RSS}{n-k-1}\) 로 정의된다. 검정통계량의 분자를 구성하는 \(TSS-RSS\) 와 분모를 구성하는 \(RSS\) 의 의미를 설명하고, 두 제곱합의 자유도를 구하는 방법도 설명하라.
개별 회귀계수의 유의성 검정에서 검정통계량을 정의하라. 귀무가설이 사실인 경우, 검정통계량의 분포는 무엇인가?
회귀모형의 결정계수는 어떻게 정의되며, 그 의미는 무엇인가?
회귀모형의 수정결정계수는 결정계수의 어떤 문제를 해결하기 위한 것인가?

2. 데이터 프레임 states는 다음과 같이 8개 변수로 구성되어 있다.

as_tibble(states) |> print(n = 3)
## # A tibble: 50 × 8
##   Population Income Illiteracy Life_Exp Murder HS_Grad Frost   Area
##        <dbl>  <dbl>      <dbl>    <dbl>  <dbl>   <dbl> <dbl>  <dbl>
## 1       3615   3624        2.1     69.0   15.1    41.3    20  50708
## 2        365   6315        1.5     69.3   11.3    66.7   152 566432
## 3       2212   4530        1.8     70.6    7.8    58.1    15 113417
## # ℹ 47 more rows

변수 Murder을 반응변수로 하고, 나머지 변수를 설명변수로 하는 회귀모형을 적합한 결과는 다음과 같다.

fit_s <- lm(Murder ~ ., states)
summary(fit_s)
## 
## Call:
## lm(formula = Murder ~ ., data = states)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4452 -1.1016 -0.0598  1.1758  3.2355 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.222e+02  1.789e+01   6.831 2.54e-08 ***
## Population   1.880e-04  6.474e-05   2.905  0.00584 ** 
## Income      -1.592e-04  5.725e-04  -0.278  0.78232    
## Illiteracy   1.373e+00  8.322e-01   1.650  0.10641    
## Life_Exp    -1.655e+00  2.562e-01  -6.459 8.68e-08 ***
## HS_Grad      3.234e-02  5.725e-02   0.565  0.57519    
## Frost       -1.288e-02  7.392e-03  -1.743  0.08867 .  
## Area         5.967e-06  3.801e-06   1.570  0.12391    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.746 on 42 degrees of freedom
## Multiple R-squared:  0.8083, Adjusted R-squared:  0.7763 
## F-statistic: 25.29 on 7 and 42 DF,  p-value: 3.872e-13

모형 fit_s에 대한 회귀모형 유의성 검정을 실시해 보자. 주어진 결과물에서 검정통계량의 값과 p-value를 구하고, 검정 결과의 의미를 해석하라.
모형 fit_s의 개별 회귀계수 검정 결과 중 5% 유의수준에서 유의적인 변수의 검정통계량 값과 p-value를 주어진 결과물에서 각각 구하라.
모형 fit_s의 결정계수와 수정결정계수를 주어진 결과물에서 각각 구하고, 그 의미를 설명하라.

3. 데이터 프레임 car_df는 230대 차량의 연비(mpg)와 배기량(displ), 변속기 종류(trans)가 입력된 자료이다.

car_df |> print(n = 3)
## # A tibble: 234 × 3
##     mpg trans  displ
##   <int> <chr>  <dbl>
## 1    18 auto     1.8
## 2    21 manual   1.8
## 3    20 manual   2  
## # ℹ 231 more rows

변수 mpg를 반응변수로 하고, 설명변수로는 연속형 변수인 displ과 범주형 변수인 trans를 사용하는 회귀모형을 적합하고자 한다. 변수 trans는 auto 또는 manual을 값으로 갖고 있다.

첫 번째 회귀모형 car_fit_1의 적합 결과는 다음과 같다. 모형 car_fit_1의 모든 회귀계수 추정값의 의미를 해석하라.

car_fit_1 <- lm(mpg ~ displ + trans, car_df)
summary(car_fit_1)
## 
## Call:
## lm(formula = mpg ~ displ + trans, data = car_df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.039 -1.357 -0.197  1.128 13.604 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  25.4609     0.5396  47.181   <2e-16 ***
## displ        -2.5520     0.1344 -18.991   <2e-16 ***
## transmanual   0.7842     0.3687   2.127   0.0345 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.548 on 231 degrees of freedom
## Multiple R-squared:  0.6446, Adjusted R-squared:  0.6415 
## F-statistic: 209.5 on 2 and 231 DF,  p-value: < 2.2e-16

두 번째 회귀모형 car_fit_2의 적합 결과는 다음과 같다. 모형 car_fit_2의 모든 회귀계수 추정값의 의미를 해석하라.

car_fit_2 <- lm(mpg ~ displ*trans, car_df)
summary(car_fit_2)
## 
## Call:
## lm(formula = mpg ~ displ * trans, data = car_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0175 -1.2676 -0.2636  1.0867 13.4728 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        25.2672     0.6311  40.036   <2e-16 ***
## displ              -2.4999     0.1606 -15.567   <2e-16 ***
## transmanual         1.3420     1.0089   1.330    0.185    
## displ:transmanual  -0.1748     0.2943  -0.594    0.553    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.552 on 230 degrees of freedom
## Multiple R-squared:  0.6451, Adjusted R-squared:  0.6405 
## F-statistic: 139.4 on 3 and 230 DF,  p-value: < 2.2e-16

두 회귀모형의 차이는 무엇인지 설명하라. 또한 주어진 결과물만을 근거로 두 모형 중 한 모형을 선택해야 한다면 어떤 모형을 선택하겠는가?

4. 데이터 프레임 auto_df에는 302대 차량의 연비와 관련된 자료가 입력되어 있다.

auto_df |> print(n = 3)
## # A tibble: 302 × 7
##     mpg cylinders displacement horsepower weight acceleration  year
##   <dbl> <fct>            <dbl>      <dbl>  <dbl>        <dbl> <dbl>
## 1    18 8                  307        130   3504         12      70
## 2    15 8                  350        165   3693         11.5    70
## 3    18 8                  318        150   3436         11      70
## # ℹ 299 more rows

모형의 평가측도에 의한 변수선택 방법 중 Best subset selection에서 최적의 변수를 선택하는 방식에 대해 설명하라.
변수 mpg를 반응변수로 하는 회귀모형을 적합하고자 한다. 설명변수는 나머지 6개 변수 중 Best subset selection으로 선택하고자 한다. 다음의 결과를 근거로 최적 변수를 선택하라.

library(leaps)
fits <- regsubsets(mpg ~ ., auto_df)

summary(fits)
## Subset selection object
## Call: regsubsets.formula(mpg ~ ., auto_df)
## 6 Variables  (and intercept)
##              Forced in Forced out
## cylinders8       FALSE      FALSE
## displacement     FALSE      FALSE
## horsepower       FALSE      FALSE
## weight           FALSE      FALSE
## acceleration     FALSE      FALSE
## year             FALSE      FALSE
## 1 subsets of each size up to 6
## Selection Algorithm: exhaustive
##          cylinders8 displacement horsepower weight acceleration year
## 1  ( 1 ) " "        " "          " "        "*"    " "          " " 
## 2  ( 1 ) " "        " "          " "        "*"    " "          "*" 
## 3  ( 1 ) " "        " "          " "        "*"    "*"          "*" 
## 4  ( 1 ) " "        " "          "*"        "*"    "*"          "*" 
## 5  ( 1 ) " "        "*"          "*"        "*"    "*"          "*" 
## 6  ( 1 ) "*"        "*"          "*"        "*"    "*"          "*"

plot(fits)

plot(fits, scale = "adjr2")