본 포스팅에서는 수식을 포함하고 있습니다.
티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니
PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다.
이번 포스팅에서는 LASSO의 명칭이 탄생하게 된 논문 'Regression Shrinkage and Selction via the LASSO'을 리뷰하고 파이썬으로 구현해보고자 한다. 여기서 다루는 내용은 다음과 같다.
3. Example -Prostate Cancer Data
4. Prediction Error and Estimation of $t$
6. Algorithms for Finding LASSO Solutions
8. Application to Generalized Regression Models
10. Results on Soft Thresholding
Summary
이 논문에서는 회귀 계수 절대값의 합에 상한을 두면서 잔차 제곱합을 최소화하는 추정법을 제안한다. 회귀 계수 절대값의 합에 상한을 둔다는 것은 $l_1$ 벌점항을 쓴다는 것이다. $l_1$ 벌점항이 추가된 최소 제곱합 문제에서 이를 만족하는 해는 특정 원소가 정확히 0이 될 수 있다는 사실이 잘 알려져 있다. 즉, Sparse 한 해를 구하는 걸로 유명하다. 이는 회귀 계수의 일부분을 정확히 0으로 추정하게 되는 것이다. 따라서 변수 선택(Subset Selection)을 해준다고 볼 수 있다. 또한 여기서 제시하는 추정 방법은 능형 회귀(Ridge Regression) 추정량처럼 분산 측면에서 안정적인 추정량이 만들어진다고 한다.
1. Introduction
기존의 최소 제곱 추정량은 다음과 같은 문제점이 있다고 한다. 첫째, 최소 제곱 추정량은 불편 추정량이지만 분산이 클 때가 종종있다. 이는 예측 오차를 크게 하여 회귀 모형의 예측 성능을 안 좋게 한다. 이런 문제는 어떻게 해결하면 좋을까? 능형 회귀처럼 회귀 추정량의 편의를 발생시키는 대신 분산을 획기적으로 줄이는 추정량을 구하는 것이 하나의 방법일 수 있다. 두 번째 문제는 모형의 해석이다. 빅데이터 시대에 데이터는 변수를 굉장히 많이 포함하고 있다. 만약 이러한 변수를 다 포함시킨 회귀 모형이 있다면 이 모형을 어떻게 해석해야 할지 막막해진다. 이 문제는 변수 선택을 통해 필요한 변수만을 모형에 추가하여 해석력을 높일 수 있다.
앞서 언급한 해결책에 대해서 좀더 살펴보자. 먼저 두 번째 문제에 대한 해결책으로 변수 선택을 이야기하였다. 그렇다면 변수 선택이 만능일까? 아니다. 변수 선택은 이산 과정(Discrete Process)이라는 문제가 있다. 왜 이산 과정이냐 하면 특정 변수가 선택되거나 빠지거나 둘 중 하나이기 때문이다. 그렇다면 이것이 왜 문제인가? 데이터가 조금만 바뀌어도 이전에 포함된 유의미한 변수가 빠질 수 있기 때문이다.
첫 번째 문제에 대한 해결방법은 대표적으로 능형 회귀가 있다. 능형 회귀는 추정량의 편의가 발생하는 대신 적절하게 회귀 계수를 축소시키고 이에 따라 분산을 작게 해주어 전체적인 예측 오차를 줄이는 방법이다. 또한 능형 회귀 추정 방법에는 튜닝 파라미터가 있는데 이 파라미터에 대하여 연속적으로 회귀 계수를 축소시켜준다(이는 변수 선택이 이산 과정인 것과 대비된다). 하지만 회귀 계수를 작게 해 줄지언정 정확하게 0으로 추정하지 못하여 변수 선택이 안된다는 단점이 있다.
저자는 이 논문 제시하는 추정 방법을 특별히 LASSO(Least Absolute Shrinkage and Selection Operator)라 칭하였고 이는 회귀 계수를 축소시켜주어 분산을 작게하고 특정 회귀 계수를 0으로 추정함으로써 변수 선택의 이점을 갖는다고 한다. 즉, LASSO는 능형 회귀와 변수 선택의 장점을 동시에 누리는 것이다.
2. The LASSO
2.1 Definition
데이터 $(\text{x}^i, y_i), i=1, \ldots, N$ 이 있다고 하자. 여기서 $\text{x}^i=(x_{i1}, \ldots, x_{ip})^t$, $y_i$는 각각 설명 변수 벡터와 반응 변수이다. 이 논문에서는 다음과 같은 가정을 한다.
가정 1)
$y_i$는 $\text{x}^i$가 주어진 경우 조건부 독립이다.
가정 2)
설명 변수 $x_{ij}$는 표준화 되어있다. 따라서 $\sum_ix_{ij}/N=1$, $\sum_ix_{ij}^2/N=1$이다.
$\hat{\beta} = (\hat{\beta}_1, \ldots, \hat{\beta}_p)^t$라 할때 LASSO 추정량 $(\hat{\alpha}, \hat{\beta})$는 다음과 같이 정의한다.
$$\DeclareMathOperator*{\argmin}{argmin} (\hat{\alpha}, \hat{\beta}) = \argmin_{\alpha\in\mathbf{R}, \beta\in\mathbf{R}^p} \sum_{i=1}^N\left(y_i-\alpha-\sum_{j=1}^p\beta_jx_{ij}\right)^2 \;\;\text{subject to }\sum_{j=1}^p|\beta_j|\leq t\tag{1}$$
여기서 $t\geq 0$ 는 튜닝 파라미터이다. 이때 모든 $t\geq 0$에 대하여 $\hat{\alpha}=\bar{y}$인 것을 알 수 있다($\alpha$는 제약식에 포함되지 않으므로 미분을 통하여 바로 구할 수 있다). 또한 반응 변수에 Centering 변환(평균을 빼주는 것)을 적용한 것을 $y_i^*$($=y_i-\bar{y}$), $\alpha^*=\alpha-\bar{y}$라 하면 (1)은 다음과 동치이다.
$$\DeclareMathOperator*{\argmin}{argmin} (\hat{\alpha}^*, \hat{\beta}) = \argmin_{\alpha^*\in\mathbf{R}, \beta\in\mathbf{R}^p} \sum_{i=1}^N\left(y_i^*-\alpha^*-\sum_{j=1}^p\beta_jx_{ij}\right)^2 \;\;\text{subject to }\sum_{j=1}^p|\beta_j|\leq t$$
따라서 일반성을 잃지 않고 $\bar{y}=0$이라 할 수 있고 이에 따라 $\alpha$를 (1)에서 제외시킬 수 있다.
$t$는 회귀 계수 $\beta_j$를 얼마나 축소시킬지 조절하는 모수이다. $\hat{\beta}^o$를 최소 제곱 추정량, $t_0=\sum_j|\hat{\beta_j^o}|$라 하자. 만약 $t<t_0$라면 (1)의 해는 최소 제곱 추정량을 축소시킨 해가 나오게 되며 $t\geq t_0$라면 (1)의 해는 최소 제곱 추정량과 같아진다. 논문에서는 $t=t_0/2$인 경우 변수의 개수가 $p/2$인 최적 모형을 찾는 것과 대략적으로 비슷해진다고 한다(뭔가 억지로 껴 맞추는 말인 거 같은데 훌륭한 분들이 썼으니 내가 모르는 의미가 있을 것 같다).
다음으로 LASSO의 모티브가 나온다. Breiman은 1993년 논문에 다음의 식을 최소화하는 Nonnegative Garrote(NG) 추정량을 제안했다.
$$\sum_{i=1}^N\left(y_i-\alpha-\sum_{j=1}^pc_j\hat{\beta}_j^ox_{ij}\right)^2 \;\;\text{subject to }c_j\geq 0, \;\sum_{j=1}^pc_j\leq t$$
Breiman은 많은 시뮬레이션을 통해 NG 추정량은 항상 변수 선택법보다 더 낮은 예측 오차를 가졌으며 실제 모형이 많은 회귀 계수가 아주 작은 양의 값을 갖는 경우를 제외하면 능형 회귀 모형과 비슷한 성능을 가질 수 있음을 보였다. NG 추정량은 최소 제곱 추정량을 포함하고 있다는 것에 주목하자. 최소 제곱 추정량은 다중 공선성이 존재한다면 좋지 않은 추정량이며 이는 NG 추정량에도 안좋은 영향을 미칠 수 있는 단점을 갖고 있다.
하지만 LASSO 추정량은 최소 제곱 추정량을 직접적으로 사용하지 않기 때문에 NG 추정량이 갖고 있는 단점을 피할 수 있다.
2.2 Orthonormal Design Case
모형 행렬(Model Matrix)를 $X=(x_{ij})$라 할 때, $X^tX=I$라고 가정하자(즉, $X$의 칼럼들은 Othonormal 이다). 이때 LASSO 추정량 $\hat{\beta}^{\text{LASSO}}$((1)의 해)는 다음과 같이 구할 수 있다.
$$\hat{\beta}^{\text{LASSO}}_j = \text{sign}\left(\hat{\beta}_j^o\right)\left(|\hat{\beta}_j^o|-\gamma\right)^+,\;\;j=1, 2, \ldots, p\tag{2}$$
여기서 $\text{sign}(x)$는 $x$가 양수이면 1 음수이면 -1 그리고 0인 경우에는 0인 함수이다. 그리고 $(x)^+=\max (0, x)$ 이다.
(2)의 증명
수식을 보면 토가 나오는 분들은 그냥 넘어가도 좋다.
(1)에서 Lagrange Multiplier를 이용하면 LASSO 추정량은 다음 식을 최소화하는 해가 된다.
$$ \frac{1}{2}\left\|X\beta-y\right\|_2^2+\gamma \|\beta\|_1\tag{2-1}$$
$\beta$를 포함하는 항을 제외하면 (2-1)을 최소화하는 문제는 아래의 식을 최소화하는 문제와 같아진다.
$$-\beta^tX^ty+\frac{1}{2}\beta^t\beta+\gamma\|\beta\|_1 \\ \\= \sum_{j=1}^p\left[\left(-\sum_{i=1}^ny_ix_{ij}\right)\beta_j+\frac{1}{2}\beta_j^2+\gamma |\beta_j|\right]\tag{2-2}$$
$X^tX=I$인 경우 최소 제곱 추정량 $\hat{\beta}^o=(X^tX)^{-1}X^ty=X^ty$이다. 따라서 $\hat{\beta}_j^o=\sum_{i=1}^ny_ix_{ij}$ 이고 (2-2)는 다음과 같이 쓸 수 있다.
$$\sum_{j=1}^p\left(-\hat{\beta}_j^o\beta_j+\frac{1}{2}\beta_j^2+\gamma |\beta_j|\right)\tag{2-3}$$
(2-3) 은 $\beta$에 대한 Separable 함수이므로 위 식을 최소화하는 $\hat{\beta}_j$를 구하면 $\hat{\beta}^{\text{LASSO}}=(\hat{\beta}_1, \ldots, \hat{\beta}_p)$ 가 된다.
(2-3)에서 $\beta_j$에 해당하는 부분을 $f(\beta_j)$라 하자. 즉,
$$f(\beta_j)=-\hat{\beta}_j^o\beta_j+\frac{1}{2}\beta_j^2+\gamma |\beta_j|$$
이다. $f(\beta_j)$는 $\beta_j$에 대한 이차함수다. 이를 최소화하는 $\hat{\beta}_j$는 $\hat{\beta}_j^o>0$인 경우와 $\hat{\beta}_j^o\leq 0$인 경우로 나누어 구할 수 있다.
case 1. $\hat{\beta}_j^o>0$
$$\hat{\beta}_j=(\hat{\beta}_j^o-\gamma)^+=\text{sign}\left(\hat{\beta}_j^o\right)\left(|\hat{\beta}_j^o|-\gamma\right)^+$$
case 2. $\hat{\beta}_j^o\leq 0$
$$\hat{\beta}_j=-(-(\hat{\beta}_j^o+\gamma))^+=\text{sign}\left(\hat{\beta}_j^o\right)\left(|\hat{\beta}_j^o|-\gamma\right)^+$$
$|\hat{\beta}_{(j)}^o|$을 $j$ 번째로 절대값이 큰 최소 제곱 회귀 추정량이라고 하자. $X^tX=I$인 상황에서 사이즈가 $k$인 최적 변수 선택은 다음과 같아진다.
$$\hat{\beta}_j^oI( \hat{\beta}_{j}^o \geq|\hat{\beta}_{(k)}^o|)\tag{3}$$
즉, Orthonormal Design인 상황인 경우 최소 제곱 추정량에서 절대값이 $\hat{\beta}_{(k)}^o$ 보다 큰 회귀 추정량들을 모두 선택하겠다는 것이다.
(3)의 증명
수식을 보면 욕이 나오는 사람들은 넘어가도 좋다.
먼저 $\mathcal{I}(\beta)=\{j : \beta_j\neq 0, \;j=1,\ldots, p\}$, $\beta^k\; (k\leq p)$를 $k$개의 0이 아닌 원소를 갖는 벡터라고 하자. 사이즈 $k\leq p$가 주어졌을 때 최적 변수 선택은 다음과 같이 잔차 제곱합을 최소화하는 $\beta^k$들을 찾는 것이다.
$$(y-X\beta^k)^t(y-X\beta^k)\tag{3-1}$$
이제 (3-1)을 정리해보자.
$$\begin{align} (y-X\beta^k)^t(y-X\beta^k) &= (\beta^k)^t\beta^k-2\beta^kX^ty+\text {const. not depend on }\beta^k \\ &= (\beta^k)^t\beta^k-2\beta^k\hat{\beta}^o+\text {const. not depend on }\beta^k \\ &=\sum_{j=1}^p\left[(\beta_j^k)^2-2\hat{\beta}_j^o\beta_j^k\right] +\text {const. not depend on }\beta^k \\ &= \sum_{j\in\mathcal{I}(\beta^k)} \left[(\beta_j^k)^2-2\hat{\beta}_j^o\beta_j^k\right]+\text {const. not depend on }\beta^k \\ &= \underbrace{\sum_{j\in \mathcal{I}(\beta^k)}(\beta_j^k-\hat{\beta}_j^o)^2+\sum_{j\in\mathcal{I}(\beta^k)}(\hat{\beta}_j^o)^2}_{\text{(3-3)}}+\text {const. not depend on }\beta^k\end{align}$$
(3-1)을 최소화하는 문제는 (3-3)을 최소화하는 문제와 같아진다. 또한 식 (3-3)은 다음의 부등식을 만족한다.
$$\text{(3-3)} \geq -\sum_{j\in\mathcal{I}(\beta^k)}|\hat{\beta}_j^o|^2 \geq -\sum_{j\in\mathcal{I}(\beta^k)}|\hat{\beta}_{(j)}^o|^2$$
첫 번째 부등식은 $\beta_j^k=\hat{\beta}_j^o$일 때 최소값을 갖는 사실을 이용하였고 두 번째 부등식은 최소 제곱 회귀 추정량의 절대값을 내림차순으로 정렬했을 경우 $k$번째까지의 합이 최소가 된다는 사실을 이용하였다. 따라서 (3)이 증명되었다.
능형 회귀 추정량은 다음 식을 최소화한다.
$$\sum_{i=1}^N\left(y_i-\sum_{j=1}^p\beta_jx_{ij}\right)+\gamma\sum{j=1}^p\beta_j^2$$
위 식을 최소화하는 능형 회귀 추정량 $\hat{\beta}^{\text{Ridge}}$은 다음과 같다.
$$\hat{\beta}^{\text{Ridge}}=(X^tX+\gamma I)^{-1}X^ty=\frac{1}{1+\gamma}\hat{\beta}_j^o, \; \gamma \geq 0\tag{4}$$
능형 회귀와 관련된 내용은 여기를 참고하자.
다음으로 NG 추정량 $\hat{\beta}_j^{\text{NG}}$은 다음과 같다.
$$\hat{\beta}_j^{\text{NG}} = \left(1-\frac{\gamma}{\hat{\beta}_j^o}\right)^+\hat{\beta}_j^o, \;\;j=1, \ldots , p\tag{5}$$
(5)의 증명
수식만 보면 분노 조절 능력이 떨어지는 사람들은 넘어가도 좋다.
NG 추정량은 먼저 다음의 식을 최소화하는 $\hat{c} = (\hat{c}_1, \ldots, \hat{c}_p)^t$를 찾는다.
$$\sum_{i=1}^N\left(y_i-\alpha-\sum_{j=1}^pc_j\hat{\beta}_j^ox_{ij}\right)^2 \;\;\text{subject to }c_j\geq 0, \;\sum_{j=1}^pc_j\leq t {\tag{5-1}}$$
NG 추정량 $\hat{\beta}_j^{\text{NG}}$은 다음과 같이 정의된다.
$$\hat{\beta}_j^{\text{NG}}=\hat{c}_j\hat{\beta}_j^o$$
논문에 나온 결과에 맞추기 위하여 식 (5-1)을 다음의 Lagrange Form으로 바꿔준다.
$$\frac{1}{2}\sum_{i=1}^N\left(y_i-\alpha-\sum_{j=1}^pc_j\hat{\beta}_j^ox_{ij}\right)^2+\gamma\sum_{j=1}^pc_j \\ =\frac{1}{2}(y-XBc)^t(y-XBc)+\gamma\sum_{j=1}^pc_j , \; \gamma \geq 0, c_j \geq 0 \tag{5-2}$$
여기서 $B$는 $j$ 번째 대각 원소가 $\hat{\beta}_j^o$인 대각 행렬이다. (5-2)에서 $c$와 상관없는 부분을 빼고 정리하면 다음과 같다.
$$\sum_{j=1}^p\left[\frac{1}{2}(\hat{\beta}_j^o)^2c_j^2-(\hat{\beta}_j^o)^2c_j+\gamma c_j\right]\tag{5-3}$$
이제 (5-3)을 최소화하는 $\hat{c}$를 찾으면 된다. (5-3)은 $c$에 대하여 Separable 함수이므로 $c$의 원소 $c_j$에 해당하는 부분만 최소화하는 문제를 생각하면 된다. 즉,
$$f(c_j) = \frac{1}{2}(\hat{\beta}_j^o)^2c_j^2-(\hat{\beta}_j^o)^2c_j+\gamma c_j\tag{5-4}$$
라고 하면 $f$를 최소화하는 $\hat{c}_j$을 찾으면 $\hat{c} = (\hat{c}_1, \ldots, \hat{c}_p)$가 (5-3)을 최소화하는 즉, 우리가 찾는 해가 된다. $f$를 정리해보자.
$$f(c_j) = \frac{1}{2}(\hat{\beta}_j^o)^2\left[c_j-\left(1-\frac{\gamma}{(\hat{\beta}_j^o)^2}\right)\right]^2 + \text{constant not depend }c_j $$
$c_j\geq 0$이기 때문에 $1-\frac{\gamma}{(\hat{\beta}_j^o)^2}>0$인 경우와 $1-\frac{\gamma}{(\hat{\beta}_j^o)^2}\leq 0$인 경우로 나누어야 한다.
case 1. $1-\frac{\gamma}{(\hat{\beta}_j^o)^2}>0$
$$\hat{c}_j = 1-\frac{\gamma}{(\hat{\beta}_j^o)^2}$$
case 2. $1-\frac{\gamma}{(\hat{\beta}_j^o)^2}\leq 0$
$$\hat{c}_j = 0$$
case 1. 과 case 2.를 종합하면 $\hat{c}_j = \left(1-\frac{\gamma}{\hat{\beta}_j^o}\right)^+$가 됨을 알 수 있다.
Orthonormal Design에서는 (2)~(5) 추정량의 특징을 비교할 수 있다. (2)~(5)를 그래프로 그려보면 다음과 같다.
(b)를 보면 능형 회귀 추정량은 최소 제곱 추정량을 축소시켜주는 것을 알 수 있다. (c)를 보면 LASSO는 최소 제곱 회귀 추정량을 축소시켜주지만 능형 회귀 추정량에서 축소시키는 방법과 다르다. 능형 회귀는 기울기를 줄이는 방식으로 축소시키지만 LASSO는 최소 제곱 추정량에서 양수를 빼주는 방식으로 축소시켜준다. 또한 LASSO는 특정 지점에서 0이 된다는 것을 알 수 있다. (d)를 보면 NG는 LASSO와 비슷해 보이지만 절대값이 큰 회귀 계수에서 LASSO보다 덜 축소된다.
2.3 Geometry of LASSO
다음 그림은 능형 회귀 추정과 LASSO 추정 방법을 기하학적으로 나타낸 것이다.
위 그림을 통해 LASSO가 왜 변수 선택이 가능한지 알 수 있다. 최소 제곱 추정량을 중심으로 하는 타원의 등고선이 사각형 영역에 닿는 점 LASSO 추정량이 되며 닿는 점이 사각형의 모서리라면 특정 회귀 계수가 0이 되므로 변수 선택이 되는 것이다. 또한 능형 회귀가 최소 제곱 추정량을 축소시켜주지만 왜 변수 선택이 안되는지 알 수 있다.
그림을 유심히 보면 다음과 같은 질문을 할 수 있다.
LASSO 회귀 추정량이 0이 아니라면 LASSO 부호와 최소 회귀 제곱 추정량의 부호가 같은가?
$p=2$이고 변수를 표준화한 경우에는 LASSO 회귀 추정량과 최소 제곱 회귀 추정량은 같은 사분면에 존재하게 된다. 따라서 LASSO 추정량이 0이 아닌 경우 최소 제곱 회귀 추정량과 부호가 같아지게 된다. 하지만 $p>2$이고 설명 변수간 상관관계가 존재한다면 상황은 달라진다. 다음 그림이 이를 보여주고 있다.
(b)를 보면 최소 제곱 추정량이 속한 곳과 LASSO 추정량이 속한 곳이 다름을 알 수 있다.
2.4 More on Two-predictor Case
$p=2$이고 $X^tX=I$라고 하자(본 논문에서는 $X^tX=I$를 언급하지 않아서 골머리 좀 썩었다). 그리고 최소 제곱 추정량 $\hat{\beta}_j^o$가 모두 양수라고 하자. 그렇다면 (2)에 의하여 LASSO 추정량이 다음과 같음을 보일 수 있다.
$$\hat{\beta}_j^{\text{LASSO}}=\left(\hat{\beta}_j^o-\gamma\right)^+$$
저자는 변수 간의 상관관계가 존재할 때 능형 회귀 추정량과 LASSO 추정량을 비교하였다. 모의실험 결과 LASSO는 상관계수의 관계없이 항상 똑같은 Solution Curve를 보여주었다고 한다. 하지만 능형 회귀 추정량은 상관계수에 따라서 다른 Solution Curve를 보인다고 한다. 또한 능형 회귀 추정량의 경우 변수간의 상관계수가 높을 때 제약식 $\sum_j\beta_j^2\leq t$ 에서 상한 $t$가 줄어듦에 따라 오히려 추정량이 축소되지 않고 오히려 커지는 경우가 발생했다고 한다. 즉, LASSO는 능형 회귀 추정 방법보다 변수간 상관관계에 대하여 Robust 하다는 것이다.
2.5 Standard Errors
LASSO 추정량은 반응 변수 $y$에 대하여 선형이 아니고 미분이 안되기 때문에 표준오차를 계산하기 어렵다. 이런 경우 붓스트랩 샘플링을 통하여 LASSO 추정량의 표준오차를 추정할 수 있다. 또한 Ridge Aprroximation을 이용하여 LASSO 추정량의 표준오차를 근사하는 방법이 있다. 이는 $|x|$가 주어진 점 $x_0$ 근방에서 다음과 같은 이차 근사식을 갖는다는 사실에 기반한다.
$$|x| \approx \frac{1}{2|x_0|}x^2+\frac{1}{2}|x_0|$$
LASSO 추정량을 특별한 언급이 없는 한 $\hat{\beta}$ 라고 하자. 위 근사식을 이용하면 다음과 같은 사실을 알 수 있다.
$$ \left\|X\beta-y\right\|_2^2+\lambda \|\beta\|_1 \approx \underbrace{\left\|X\beta-y\right\|_2^2+\frac{1}{2}\lambda\sum_{j=1}^p\frac{\beta_j^2}{|\hat{\beta}_j|}}_{\text{(6-1)}}+c\tag{6}$$
여기서 $c$는 $\beta_j$와 관계없는 상수이다. 따라서 근사식을 이용하면 (6)의 왼쪽을 최소화하는 문제는 (6-1)을 최소화하는 문제로 바꿀 수 있다. (6-1)은 $\beta$에 대한 이차식이며 미분이 가능하므로 Analytic한 해를 찾을 수 있다. 이를 $\hat{\beta}^*$라 하면 다음과 같음을 알 수 있다.
$$\hat{\beta} \approx \hat{\beta}^* = (X^tX+\lambda W)^{-1}X^ty\tag{7}$$
여기서 $W$은 $j$ 번째 대각 원소가 $1/|\hat{\beta}_j|$인 대각 행렬이다. 이제 (7)을 이용하면 LASSO 추정량의 표준오차를 다음의 분산 추정량을 이용하여 추정할 수 있다.
$$\widehat{\text{var}}(\hat{\beta})\approx\widehat{\text{var}}(\hat{\beta}^*)=\hat{\sigma}^2(X^tX+\lambda W)^{-1}X^tX(X^tX+\lambda W)^{-1}\tag{8}$$
여기서 $\hat{\sigma}^2$은 $y$의 분산 추정량이다. 하지만 $\hat{\beta}_j=0$인 경우 $\hat{\beta}_j$에 대한 분산을 0으로 추정한다는 단점이 있다고 한다. Ridge Approximation을 통해 Iterative Ridge Estimation 방법을 쓰면 꽤 괜찮은 LASSO 추정량을 얻을 수 있을 것 같지만 굉장히 비효율적이라고 한다. 그래도 제약식의 상한 $t$를 추정하는데 이 방법이 유용하다고 한다.
3. Example - Prostate Cancer Data
논문에서는 Prostate Cancer Data를 이용하여 최소 제곱 모형, LASSO 모형, 그리고 변수 선택법을 이용한 모형을 적합하여 각 추정량을 비교하였다. 아래 그림은 $s=t/\sum_{j=1}^p |\hat{\beta}_j^o|$의 변화에 따른 LASSO 회귀 추정량을 그래프로 나타낸 것이다.
위 그림에서는 LASSO 회귀 추정량이 $s$가 0에 다가갈수록 회귀 추정량이 음수인 경우에는 단조롭게 증가하고 양수인 경우에는 단조롭게 감소한다. 하지만 이러한 단조성은 항상 만족하지 않는다. 이는 아래 그림을 보면 직관적으로 알 수 있다.
위 그림을 보면 주어진 $s$에 대하여 최소 제곱 추정량 $\hat{b}_2$과 LASSO 추정량 $b_2(s)$이 모두 양수이지만 $\hat{b}_2<b_2(s)$ 인 것을 알 수 있다. 따라서 이런 경우에는 단조성이 만족하지 않을 수 있다.
아래 Table 1은 모형 적합 결과를 나타낸 것이다.
먼저 변수 선택과 최소 제곱 모형을 비교해보자. 선택된 변수에 대하여 회귀 추정치, Z-score가 변수 선택 모형에서 더 높다(절편항 제외). 이는 변수간 양의 상관관계가 있는 경우 종종 일어나는 일이라고 한다. 그리고 LASSO 결과는 $\hat{s}=0.44$에서 계산된 회귀 추정치이다. $\hat{s}$는 교차 검증을 통하여 구했다고 한다. 이에 대한 내용은 뒤에 나온다. LASSO와 최소 제곱 모형을 비교해보면 절편항은 같고 나머지 회귀 추정치가 LASSO에서 더 작다는 것을 알 수 있다. 즉, 축소가 잘 된 것이다. 이러한 축소 현상이 Z-score 에도 영향을 미쳐서 LASSO Z-score가 더 작아졌다. 사실 Z-score가 더 작아지는 것은 직관적으로 잘 모르겠다. 왜냐하면 축소된 추정량은 분산도 작아지기 때문이다(저자는 이부분에 대한 언급은 없다).
맨 끝에서 2번째 칼럼에는 LASSO 회귀 추정량의 분산 추정치가 정리되어있다. 이는 $\hat{s}=0.44$인 경우에 붓스트랩 샘플링을 이용한 표준 오차 추정량이다. 아래 표는 붓스트랩 샘플링과 Ridge Approximation (8)을 이용한 표준 오차 근사 추정치를 계산하여 정리한 것이다.
위 표를 보면 0이 아닌 회귀 추정치에 대하여 Ridge Approximation을 이용한 표준 오차가 고정된(Fixed) $t$에 대한 붓스트랩 샘플링 추정량을 잘 근사하는 것으로 보인다. 아래 그림은 $\hat{s}=0.44$에 대하여 200개의 붓스트랩 샘플을 이용하여 회귀 추정치를 구한 것이다.
붓스트랩 샘플을 이용한 90% 신뢰 구간을 구해본 결과 lcavol, svi를 제외하고는 모두 0을 포함했다고 한다. LASSO 회귀 추정치가 0인 회귀 계수의 신뢰구간은 모두 0을 포함했다. 즉, 의미 없다고 판단된 변수는 신뢰구간을 통해서도 의미 없음이 증명된 것이다. 또한 이 결과는 데이터가 조금 바뀌더라도 변수 선택 결과는 거의 일치한다는 것을 암시한다.
4. Prediction Error and Estimation of $t$
여기서는 제약식의 상한 $t$를 선택하는 방법에 대하여 알아본다. 논문에서는 교차 검증, 일반 교차 검증, 그리고 Risk의 불편 추정을 통하여 Analytical 하게 구할 수 있는 방법을 소개한다. 첫 번째, 두 번째 방법은 $X$가 랜덤인 경우에 적용할 수 있는 방법이며 세 번째 방법은 $X$가 고정된 상수일 때 적용 가능한 방법이다. 저자는 실제 사례에서는 $X$가 랜덤인 경우와 고정된 상수인 경우가 명확히 구분(?)되지 않으므로 그냥 편한 거 쓰면 된다고 한다. Experimental Study와 Observational Study에서 전자는 주로 X가 고정이고 후자는 랜덤으로 본다. 그리고 그렇게 해야 하는 근거는 꽤나 명확하다. 근데 왜 명확한 구분이 없다고 했는지 모르겠다. 훌륭하신 분들이 그렇다니까 그냥 그런가 보다 하고 넘어가야지...
$X$와 $Y$의 관계가 다음과 같다고 가정하자.
$$Y = \eta (X)+\epsilon$$
여기서 $E(\epsilon)=0$, $\text{var}(\epsilon)=\sigma^2$이다. 우리가 여러 모형 적합 방법을 이용하여 추정한 $\eta (X)$를 $\hat{\eta}(X)$라 할 때 $\hat{\eta}(X)$의 평균 제곱 오차(Mean Square Error : ME)는 다음과 같이 정의한다.
$$\text{ME} = E\left\{\hat{\eta}(X)-\eta (X)\right\}^2$$
그리고 예측 오차(Prediction Error : PE)는 다음과 같이 정의한다.
$$\text{PE} = E\left\{Y-\hat{\eta}(X)\right\}^2 = \text{ME} + \sigma^2 \tag{9}$$
교차 검증은 예측 오차를 추정하는 방법이다. 이때 $s=t/\sum\hat{\beta}_j^o$는 0과 1 사이의 값을 쪼개서 교차 검증(예 : 5-fold 교차 검증)을 진행하고 가장 작은 예측 오차를 추정하는 $s$를 최적 $\hat{s}$로 선택한다.
다음으로 일반 교차 검증을 이용하여 $t$값을 추정하는 방법을 소개한다.
벌점항 $\lambda\sum|\beta_j|$를 $\lambda\sum\beta_j^2/|\beta_j|$로 바꿔서 생각한다면 바꿔진 벌점항에 대한 추정량 $\tilde{\beta}$는 다음과 같이 구할 수 있을 것이다.
$$\tilde{\beta}=(X^tX+\lambda W^-)^{-1}X^ty\tag{10}$$
여기서 $W=\text{diag}(|\hat{\beta}_j|)$, $W^-$은 $W$의 일반화 역행렬(Generalized Inverse)이다. 이때 효율적인 파라미터 개수(the number of effective parameter)는 다음과 같이 근사 할 수 있다.
$$p(t) = \text{tr}\{X(X^tX+\lambda W^-)^{-1}X^t\}$$
$\text{rrs}(t)=(X\tilde{\beta}-y)^t(X\tilde{\beta}-y)$라 하자. 즉, $\text{rrs}(t)$는 잔차제곱합이다. 이때 다음과 같이 일반 교차 검증 통계량 $\text{GCV}(t)$를 정의한다.
$$\text{GCV}(t) = \frac{1}{N}\frac{\text{rrs}(t)}{\{1-p(t)/N\}^2}\tag{11}$$
이때 $\text{GCV}(t)$를 최소화하는 $t$를 선택한다.
마지막으로 추정량의 Risk를 이용하여 최적 $t$를 선택하는 방법을 소개하고 있다. $g : \mathbf{R}^p\rightarrow\mathbf{R}^p$인 almost differential function이라고 할 때 $\mu$의 추정량 $\hat{\mu}$은 다음과 같이 표현했다고 가정하자.
$$\hat{\mu} = z + g(z)$$
여기서 $z\sim N(\mu, I)$인 확률변수이다. 이때 Stein (1988)은 다음과 같은 사실을 증명하였다.
$$E\|\hat{\mu}-\mu\|^2 = p + E\left(\|g(z)\|^2+2\sum_{j=1}^pdg_j/dz_j\right)\tag{12}$$
저자는 $\hat{\beta}_j^o/\hat{\tau}$이 근사적으로 표준 정규분포를 따른다는 사실과 (12)를 이용하여 다음과 같이 평균 제곱 오차의 불편 추정량을 근사적으로 제시하였다.
$$R\{\hat{\beta}(\gamma)\} \approx \hat{\tau}^2\left\{ p-2\cdot \text{#} (j : |\hat{\beta}_j^o/\hat{\tau}|\leq\gamma ) +\sum_{j=1}^p \min (|\hat{\beta}_j^o/\hat{\tau}|,\gamma)^2 \right\}\tag{13}$$
논문에는 min이 아닌 max라고 나왔는데 틀린 표현이다.
(13) 증명
먼저 $\hat{\mu}=(\hat{\mu}_1,\ldots ,\hat{\mu}_p), \;\hat{\mu}_j = \text{sign}(x_j)(|x_j|-\gamma)^+, \;\; x_j\sim N(\mu_j, 1)$인 경우 $E\|\hat{\mu}-\mu\|^2$의 추정량은
$$p-2\cdot \text{#} (j : |x_j|\leq\gamma ) +\sum_{j=1}^p \min (|x_j|,\gamma)^2\tag{13-1}$$
임을 먼저 보일 것이다.
$$\begin{align} \hat{\mu}_j &= \text{sign}(x_j)(|x_j|-\gamma)^+ \\ \\ &= x_j + \underbrace{\text{sign}(x_j)(|x_j|-\gamma)^+ - x_j}_{g_j(x_j)}\end{align}$$
여기서 다음과 같은 사실을 알 수 있다.
$$\begin{align}g_j(x_j)^2 &= \gamma^2 & \text{if } |x_j|>\gamma \\ \\ &= x_j^2 & \text{o.w} \end{align}$$
그러므로 $g_j(x_j)^2 = \min (|x_j|, \gamma)^2$이다. 또한
$$\begin{align}dg_j/dx_j &= 0 & \text{if } |x_j|>\gamma \\ \\ &= -1 & \text{o.w}\end{align}$$
이다. 따라서
$$\|g(x)\|^2+2\sum_{j=1}^pdg_j/dx_j = \sum_j^p\min (|x_j|, \gamma)^2) -2\cdot \text{#}(j : |x_j|\leq \gamma)$$
이고 (12)의 결과를 이용하면 (13-1)을 얻을 수 있다.
한편 $x_j = \hat{\beta}_j^o/\hat{\tau}$라고 하면 $x_j \approx N(0,1)$이다. $\hat{\beta}_j/\hat{\tau} = \hat{\mu}_j$라 하면 다음과 같음을 알 수 있다.
$$R(\hat{\beta}(\gamma))\approx \hat{\tau}^2(p+\sum_j^p\min (|x_j|, \gamma)^2-2\cdot \text{#}(j : |x_j|\leq \gamma))$$
저자는 $R(\hat{\beta}(\gamma))$를 가장 작게 하는 $\gamma$를 선택하라고 한다. 이때 선택한 값을 $\hat{\gamma}$이라 하면 이에 대응하는 $t$를 다음과 같이 추정한다.
$$\hat{t} = \sum_{j=1}^p(|\hat{\beta}_j^o|-\hat{\gamma})^+$$
저자는 Stein 방법을 이용한 $t$의 추정 방법이 계산상에서 훨씬 간편하다고 한다.
5. LASSO as Bayes Estimate
여기서는 LASSO 추정방법을 Bayesian 추정 관점으로 해석하였다. LASSO는 $\beta_j$에 대한 사전 분포를 double exponential 분포로 하고 Posterior mode를 구하는 것으로 볼 수 있다.
6. Algorithms for Finding LASSO Solutions
주어진 $t\geq 0$에 대하여 (1)은 $2^p$개의 부등식 제약 조건이 있는 최소 제곱 문제를 푸는 것으로 표현할 수 있다. $2^p$는 $\beta_j$의 부호의 조합을 의미한다. Lawson과 Hansen은 부등식 제약 조건 $G\beta \leq h$가 있는 상황에서 최소 제곱 문제를 푸는 방법을 제공했다. 이때 $G$는 $m\times p$ 행렬이고 $\leq$는 원소별로 부등식을 적용한 것이다. LASSO의 경우 $m=2^p$인 상황이다. 따라서 $p$가 커질수록 행렬 계산량이 엄청나게 늘어나서 실용적이진 않다.
$2^p$개의 제약 부등식을 모두 고려하는 것이 아니라 순차적으로 고려하여 Kuhn-Tucker 조건을 만족하는 feasible 한 솔루션을 찾아가는 방식으로 문제를 풀 수 있다. $g(\beta)=\sum_{i=1}^N(y_i-\sum_j\beta_jx_{ij})^2$, $\delta_i, \;i=1, 2, \ldots, 2^p$을 $(\pm 1, \ldots, \pm 1)$ 형식의 $p$차원 벡터라고 하자. 그렇다면 $\sum |\beta|_j \leq t$ 조건은 $\delta_i \leq t, \; \forall i=1, \ldots, 2^p$와 동치이다. 그렇다면 LASSO 추정량을 구하는 알고리즘은 다음과 같다.
알고리즘
(1) $E=\{i_0\}$, $\delta_{i_0} = \text{sign}(\hat{\beta}^o)$라 하자. $\hat{\beta}^o$는 최소 제곱 추정량이다.
(2) 만약 $\sum |\hat{\beta}_j^o| \leq t$ 라면 알고리즘을 종료한다.
(3) 인덱스 $\{i : \delta_i\hat{\beta}_j^o > t\}$을 하나 찾은 다음 $E$에 추가한다.
(4) 다음을 만족하는 $\hat{\beta}$를 찾는다.
$$\text{minimize } g(\beta) \;\;\text{subject to } G_E\beta \leq t\mathbf{1}$$
여기서 $G_E$는 $\delta_i \; i\in E$를 행으로 하는 행렬이고 $\mathbf{1}$은 $G_E$의 행 벡터 길이와 같은 모든 원소가 1인 벡터이다.
(5) 만약 $\sum |\hat{\beta}_j| \leq t$ 라면 알고리즘을 종료한다.
(6) 인덱스 $\{i : \delta_i\hat{\beta} > t\}$을 하나 찾은 다음 $E$에 추가하고 (4)로 돌아간다.
이 알고리즘은 유한 스텝을 거치면 반드시 종료하게 된다. 왜냐하면 한 스텝에 하나의 원소가 $E$에 들어가게 되고 이렇게 들어갈 수 있는 원소는 $2^p$로 유한하기 때문이다. 이렇게 하여 얻어진 솔루션은 우리가 찾고자 하는 LASSO 추정량이 된다. 저자는 많은 실험 결과 $(0.5p, 0.75p)$사이에 횟수를 거치면 알고리즘이 종료되는 것을 발견했다고 한다(물론 항상 그런 것은 아니다). 따라서 $p$가 커질수록 엄청난 계산과 시간이 필요하게 된다. LASSO 추정량을 찾기 위한 여러 효율적인 알고리즘이 있다 하니 관심 있는 사람들은 찾아보자.
나는 위의 알고리즘을 파이썬으로 구현하고 Prostate Cancer Data에 적용해보았다. 먼저 Prostate Cancer Data를 여기서 다운받았다.
필요한 모듈을 임포트하고 데이터를 불러오자.
import numpy as np
import pandas as pd
import statsmodels.api as sm
import time
import matplotlib.pyplot as plt
import seaborn as sns
from functools import reduce
from itertools import product
from tqdm import tqdm
from scipy.linalg import svd
from scipy.optimize import nnls
from random import choices, seed
df = pd.read_csv('prostate_cancer.csv')
논문에서 제시한 알고리즘의 핵심은 LASSO는 부등식 제약이 있는 최소 제곱 문제(Least Square with Inequality Constraints : LSIC)를 풀어야 한다는 것이다. 논문에서 인용한 문헌 Charles L. Lawson, Richard J. Hanson - Solving Least Squares Problems(1987)
을 확인해보니 LSIC를 Least Distance Problem(LDP)으로 바꿔서 풀었다. 또한 LDP는 Non Negative Least Square(NNLS) 알고리즘을 활용한다는 사실도 알았다. 이를 토대로 LASSO 추정량을 구하는 알고리즘을 만들어보았다.
먼저 LDP를 풀어주는 함수를 만들었다. 이 함수는 다음과 같은 문제를 풀어낸다.
$$\text{Minimize } \|x\| \;\; \text{subject to }Gx \geq h$$
def ldp(G,h):
E = np.vstack((G.T,h))
n = G.shape[1]
f = np.hstack((np.zeros(n), [1]))
u_hat, _ = nnls(E, f)
residual = np.matmul(E, u_hat) - f
norm_residual = np.linalg.norm(residual)
if norm_residual == 0:
print('not feasible')
return False
else:
x_hat = []
for j in range(n):
x_hat.append(-residual[j]/residual[n])
return x_hat
위 코드에서 5 번째 줄을 보면 nnls함수가 있는데 이는 NNLS를 풀어주는 함수이다. 원래 이것도 직접 구현하였으나 특정 부분에서 에러가 발생했고 이 당시 피곤해서 더 고치고 싶지 않았다. scipy에서 제공하는 nnls 함수를 사용했다. nnls 함수는 다음과 같은 문제를 풀어낸다.
$$\text{Minimize } \|Ex-f\|^2 \;\;\text{subject to }x_j\geq 0, \; \forall j$$
그리고 아래 함수가 가장 중요한 함수인 LSIC를 풀어내는 함수이다. 이 함수는 다음의 문제를 풀어낸다.
$$\text{Minimize } \|Ex-f\|^2 \;\;\text{subject to }Gx \geq h$$
def ls_inequality_constraints(E,f,G,h):
U, S, V_t = np.linalg.svd(E)
V = V_t.T
num_zeros = U.shape[0] - len(S)
num_non_zeros = len(S)
S = np.diag(S)
temp = np.matmul(U.T,f)
f1 = temp[:num_non_zeros]
S_inv = np.linalg.inv(S)
tilde_G = reduce(np.dot, [G, V, S_inv])
tilde_h = h - np.matmul(tilde_G, f1)
z = ldp(tilde_G, tilde_h)
x_hat = reduce(np.dot, [V, S_inv, z+f1])
return x_hat
다음으로 몇 가지 필요한 함수를 정의했다. 칼럼을 표준화시켜주는 함수와 벡터의 부호를 출력하는 함수이다.
def standardize_column(x):
mean = np.mean(x)
std = np.std(x)
return (x-mean)/std
def sign(x):
res = []
for v in x:
if v > 0:
res.append(1)
else:
res.append(-1)
return res
마지막으로 LASSO 추정량을 찾아주는 함수이다.
def lasso(X,y,t):
if isinstance(X, pd.DataFrame):
X = X.values
fit = sm.OLS(y,X).fit()
beta = np.array(fit.params) ## initial vector
sign_vectors = []
num_parameter = len(beta)
for i in range(num_parameter):
sign_vectors.append([1,-1])
beta_sign = tuple(sign(beta))
G = np.expand_dims(np.array(beta_sign), axis=0)
sign_set = set(product(*sign_vectors))
sign_set = sign_set-{beta_sign}
while True:
constraint_state = True
for x in list(sign_set):
temp_x = np.array(x)
if temp_x.dot(beta) > t:
constraint_state = False
G = np.vstack((G,temp_x))
sign_set = sign_set-{x}
break
if constraint_state:
return beta
h = t*np.ones(G.shape[0])
beta = ls_inequality_constraints(X,y,-G,-h)
이제 Prostate Cancer Data에 적용해보자. 먼저 설명 변수를 표준화시키고 최소 제곱 모형을 적합한다. 그런 다음 회귀 계수 절대값의 합을 구해준다.
train_x = df[['lcavol','lweight','age','lbph','svi','lcp','gleason','pgg45']]
train_y = df['lpsa']
X = train_x.apply(standardize_column)
y = train_y
fit = sm.OLS(y,X).fit()
sum_abs_coef = sum(abs(fit.params))
다음 코드는 $s$의 값을 정해둔 뒤 이를 $t=s\cdot\sum_j|\hat{\beta}_j^o|$로 바꾸어준다. 그리고 각 $t$에 대하여 LASSO 추정량을 구한다.
s = np.linspace(0,1,101)
ts = s*sum_abs_coef
beta_lasso = []
for t in tqdm(ts):
beta_lasso.append(lasso(X,y,t))
beta_lasso = np.array(beta_lasso)
LASSO 추정량을 시각화해보자.
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
for i in range(len(train_x.columns)):
plt.plot(s,beta_lasso[:,i],label=train_x.columns[i])
plt.axvline(0.44,linestyle='--',color='k',linewidth=1)
plt.legend()
plt.show()
추정량의 궤적이 조금 다르지만 거의 비슷한 거 같다 ㅎㅎ;;
그리고 최종 LASSO 추정량과 붓스트랩 샘플링을 이용하여 표준오차를 계산해보았다. 여기서는 Fixed $t\;(s=0.44)$에 대해서 구했다.
seed(1)
s_hat = 0.44
t_hat = s_hat*sum_abs_coef
beta_hat = lasso(X,y,t_hat)
for j in range(len(beta_hat)):
if beta_hat[j] < 10**(-15):
beta_hat[j] = 0
res_df = pd.DataFrame()
res_df['Predictor'] = ['Intercept']+list(train_x.columns)
res_df['Coefficient'] = np.hstack((y.mean(),beta_hat))
B = 200
bootstrap_betas = []
for i in tqdm(range(B)):
bootstrap_residual = choices(fit.resid, k=len(fit.resid))
bootstrap_residual = np.array(bootstrap_residual)
bootstrap_y = fit.fittedvalues + bootstrap_residual
bootstrap_intercept = np.mean(bootstrap_y)
bootstrap_beta = lasso(X, bootstrap_y, t_hat)
bootstrap_betas.append(np.hstack((np.mean(bootstrap_y),bootstrap_beta)))
bootstrap_betas = np.array(bootstrap_betas)
boostrap_std = []
for i in range(len(train_x.columns)+1):
std = np.std(bootstrap_betas[:,i],ddof=1)
boostrap_std.append(std)
res_df['Standard Error'] = boostrap_std
결과를 비교해보자.
변수 선택은 동일하게 되었으나 회귀 계수에 조금 차이가 난다. 그리고 표준 오차의 경우 대부분은 비슷하나 Age 변수에 대한 회귀계수의 표준오차가 논문에서와 큰 차이가 발생했다. ㅠ.ㅠ
다음으로 박스 플롯을 그려보았다. 박스 플롯을 그리는 방법은 여기를 참고하기 바란다.
## 박스 그림
num_group = len(train_x.columns) ## 그룹개수
edgecolor = 'blue' ## 선 색깔
linewidth = 2 ## 선 굵기
facecolors = sns.color_palette('hls',num_group) ## 면색
medianlinecolor = 'white' ## 중앙선 색
data = []
for i in range(1, len(train_x.columns)+1):
data.append(bootstrap_betas[:,i])
fig = plt.figure(figsize=(8,8)) ## 캔버스 생성 사이즈 8*8
fig.set_facecolor('white') ## 캔버스 배경색 설정
box = plt.boxplot(data, patch_artist=True) ## 상자 수염 그림 출력
## 상자, 1.5*IQR 초과/미만 데이터(flier), 중앙선 넘어가는 데이터 꾸미기
for item in ['boxes', 'fliers', 'medians']:
for i, v in enumerate(box[item]):
if item == 'boxes':
plt.setp(v, facecolor=facecolors[i], edgecolor=edgecolor)
elif item == 'fliers':
plt.setp(v, markerfacecolor=facecolors[i])
else:
plt.setp(v, color=medianlinecolor)
## 최상단, 최하단 선분(caps)과 상자와 이어지는 선분(whisker) 꾸미기
for item in ['whiskers','caps']:
for i, v in enumerate(zip(box[item][::2],box[item][1::2])):
if item == 'whiskers':
plt.setp(v, color=edgecolor, linewidth=linewidth, linestyle='--')
else:
plt.setp(v, color=edgecolor, linewidth=linewidth)
fontsize = 13 ## 폰트 사이즈
plt.xticks(range(1,len(train_x.columns)+1),train_x.columns, fontsize=fontsize) ## x축 눈금 라벨
plt.title('Box Plot',fontsize=fontsize) ## 타이틀
plt.show()
역시 논문과 흡사하게 나왔다. 이걸로 만족해야지 ㅋㅋㅋㅋ
7. Simulations
7.1 Outline
저자는 모의실험을 통하여 최소 제곱 추정량과 LASSO, NG, 변수 선택, 능형 회귀의 추정량들을 비교하였다.
7.2 Example 1
실제 모형이 어느 정도 Sparse 한 경우 LASSO의 성능이 가장 좋았다고 한다. LASSO가 실제 파라미터를 정확하게 추정한 경우, 즉 0이 아닌 계수는 잘 추정하고 나머지는 0으로 추정한 경우는 2.5%에 그쳤지만 실제 파라미터를 포함하는 경우는 95.5%였다고 한다. 반면 변수 선택으로 인하여 실제 파라미터를 정확하게 선택하는 24%로 LASSO 보다 더 높다. 즉, 0이 아닌 회귀 계수를 정확하게 선택하는 능력은 변수 선택법이 LASSO보다 더 좋다는 것이다.
7.3 Example 2
여기서는 모든 계수가 0.8인 경우로 다시 말하면 모든 계수가 0이 아닌 실제 모형을 이용하여 모의실험을 하였다. 이 경우에는 능형 회귀가 가장 좋은 성능을 보였다고 한다. LASSO가 그 뒤를 이었다. 즉, Sparse 한 모형이 아닌 경우 능형 회귀가 더 좋을 수 있다는 것을 보여주고 있다.
7.4 Example 3
여기서는 8개의 계수중 하나만 0이 아닌 경우에 대하여 모의실험을 하였다. 여기서는 Garotte 추정량과 변수 선택법이 가장 좋은 성능을 보였고 그 뒤를 LASSO가 따르고 있다. 즉, 굉장히 Sparse 하면 LASSO 보다 변수 선택법과 Garrote 추정량이 예측오차면에서 좋다고 할 수 있다. 또한 이 경우 능형 회귀가 가장 좋지 않은 성능을 보였는데 굉장히 Sparse하면 능형 회귀는 매우 좋지 않음을 보여주는 예이다.
7.5 Example 4
이번에는 계수가 40개인 꽤나 큰 모형에 대하여 실험을 진행하였다. 40개 중에 0이 아닌 계수는 절반인 20개였다. 이 경우 능형 회귀가 가장 좋은 성능을 보였으며 다음으로 LASSO가 좋은 성능을 보였다고 한다.
전체적으로 보면 어느 경우에서나 LASSO가 상위권에 있어서 쓸만한 추정방법이라고 얘기하는 것 같다. 하지만 항상 그러하듯이 언제나 좋은 것 또한 아닌 것을 보여준다.
8. Application to Generalized Regression Models
여기서는 LASSO의 아이디어를 일반화 모형에 적용할 수 있음을 보여주었다. 예를 들어 일반화 선형 모형의 경우 최대 우도 추정량을 구하게 되는데 우도 함수가 복잡하여 Newton Raphson 알고리즘과 같은 Iteratively Reweighted Least Square(IRLS) 방법을 이용하여 추정하게 된다. 이때 저자는 IRLS에서 각 스텝마다 LASSO를 적용하여 추정량을 구할 수 있다고 하였다. 비록 추정량의 수렴성이 일반적으로 보장되지는 않더라도 저자의 경험으로는 꽤나 잘 동작했다고 한다.
8.1 Logistic Regression
저자는 로지스틱 회귀 모형을 LASSO를 이용하여 추정하였다. 왜냐하면 로지스틱 회귀 모형은 IRLS를 통해서 추정해야 하기 때문이다. 여기서 Backward Stepwise를 이용하여 최종적으로 선택한 모형과 LASSO를 이용한 모형을 추정하였다. LASSO 추정량은 5번 정도의 스텝을 거쳐 $10^{-6}$의 범위의 수렴성을 보여주었다고 한다.
9. Some Further Extensions
현재 저자는 LASSO의 아이디어를 두 가지 분야에 적용시키는 것을 연구하고 있다고 한다. 하나는 의사 결정 나무 모델이고 다른 하나는 Multivariate Adaptive Regression Splines(MARS)라고 한다. LASSO와 같은 축소 과정은 의사 결정 나무의 가지치기보다 더 정확하고 해석력 좋은 결정 나무를 만들어낼 수 있다고 한다. 또한 저자는 MARS에 LASSO 아이디어를 접목하여 모형을 키우거나 축소시키는 알고리즘을 개발 중이라고 한다. 또한 모형 행렬이 Full Rank가 아닌 문제에 대해서도 LASSO 아이디어를 적용시킬 수 있다고 한다.
10. Results on Soft Thresholding
여기서는 Orthonormal Design인 경우의 LASSO 추정량 $\text{sign}\left(\hat{\beta}_j^o\right)\left(|\hat{\beta}_j^o|-\gamma\right)^+$가 가장 이상적인 변수 선택자의 성능면에서 밀접한 연관이 있다고 한다. 저자는 LASSO 추정량에 대한 Risk와 이상적인 변수 선택자의 Risk의 관계를 소개하였다. 이는 LASSO의 잠재적인 유용성을 암시한다고 한다. 하지만 실제 문제에서 Orthonormal Design인 경우는 거의 없고 대부분 변수간 상관관계가 존재한다. 이런 상황에서 위와 비슷한 결과를 얻어내는 것이 힘들어 보인다고 한다.
11. Discussion
이 논문에서는 연속적인 축소 연산을 이용하여 특정 회귀 계수를 정확하게 0으로 추정할 수 있음을 보였다. 이를 통해 LASSO는 변수 선택과 능형 회귀의 장점을 갖추었고 이들의 경쟁자가 될 수 있음을 보였다. 저자는 다음 세 가지 경우에 대하여 LASSO, 변수 선택, 능형 회귀 방법 중 어떤 방법이 상대적으로 좋을지를 조사하였다.
(a) 큰 영향을 미치는 적은 수의 변수 - 변수 선택법이 가장 좋은 성능을 보였으며, LASSO는 별로 좋지 않고 능형 회귀는 가장 안 좋다고 한다.
(b) 적당한 영향을 미치는 중간 정도 숫자의 변수 - LASSO가 가장 좋은 성능을 보였으며 능형 회귀, 변수 선택 순으로 좋았다고 한다.
(c) 약한 영향을 미치는 많은 숫자의 변수 - 능형 회귀, LASSO 그리고 변수 선택 순으로 좋았다고 한다.
Garotte 추정량은 (a)에서는 LASSO 보다 조금 좋았고 (b)에서는 LASSO 보다 조금 안 좋았다고 한다. 저자는 성능적인 측면을 제외하더라도 변수 선택, LASSO, Garotte 추정방법이 능형 회귀보다 더 해석이 좋은 모형을 만들어낸다는 장점이 있다고 한다.
댓글