[회귀 분석] 10. 가중 최소 제곱법(Weigted Least Square)으로 회귀 모형 적합하기 with Python

---

안녕하세요~ 꽁냥이에요!

 

선형 회귀 모형의 가정 중에서 오차가 설명변수에 의존하지 않는 등분산성 가정이 있습니다. 하지만 때때로 이 가정을 만족하지 않는 상황이 발생할 수 있는데요. 이런 상황에서 최소 제곱 회귀 추정량은 좋지 않은 성질을 갖고 있지요.

 

따라서 최소 제곱법이 아닌 다른 추정법을 이용하여 회귀 모형을 만들어야 합니다. 이때 사용하는 것이 오늘 소개할 가중 최소 제곱법(Weighted Least Square)입니다. 

 

이번 포스팅에서는 가중 최소 제곱법(Weighted Least Square)과 파이썬(Python)을 이용한 예제를 알아보겠습니다.

 

1. 가중 최소 제곱법(Weighted Least Square)이란?

2. 가중 최소 제곱(Weighted Least Square) 모형 적합하기 with Python

3. statsmodels 모듈을 이용하여 가중 최소 제곱 모형 적합하기

---

   1. 가중 최소 제곱법(Weighted Least Square)이란?

- 정의 -

먼저 우리에게 데이터 $(\tilde{x}_i, y_i), \;(i=1, \cdots, n)$ 가 있다고 합시다. 여기서 $\tilde{x}_i=(x_{i1}, x_{i2}, \cdots, x_{ip})^t$ 입니다. 이때 실제(true) 선형 회귀 모형이 다음과 같다고 할게요.

$$y_i = \beta_0 + \beta_1x_{i1} + \cdots + \beta_px_{ip} + \epsilon_i, \:\: i=1, \ldots, n$$

 

이 때 다음과 같이 가정합니다.

---

가정 1. $x_{i1}, x_{i2}, \ldots, x_{ip}$는 상수(Fixed Constant)라고 가정합니다.

가정 2. $\epsilon_i$는 서로 독립이며 정규분포 $N(0, \sigma_i^2)$을 따릅니다.

---

(가정 2를 살펴보면 오차의 등분산성을 만족하지 않는 상황도 포함되어 있다는 것을 알 수 있습니다.)

 

$\beta=(\beta_0, \cdots, \beta_p)^t$라고 할 때 가중 최소 제곱법은 다음의 가중 오차 제곱합 $L(\beta)$을 최소화시키는 회귀모형을 찾는 방법입니다.

$$L(\beta)=\sum_{i=1}^nw_i(y_i-\beta_0-\beta_1x_{i1}-\cdots-\beta_px_{ip})^2\tag{1}$$

 

여기서 $w_i >=0\;(i=1, \ldots, n)$입니다. 식 (1)을 행렬로 표현하면 아래와 같습니다.

 

$$L(\beta) = (y-X\beta)^tW(y-X\beta)$$

 

여기서 $y=(y_1, \cdots, y_n)^t$, $W$는 $w_i$를 대각 원소로 하는 $n\times n$ 대각 행렬, 그리고 $X$는 상수항과 설명변수로 이루어진 $n\times (p+1)$ 행렬입니다.

 

$$\begin{equation}
X=
\begin{pmatrix} 1 & x_{1,1} & \cdots & x_{1,p} \\ 1& x_{2,1} & \cdots & x_{2,p} \\ \vdots & \vdots & \ddots & \vdots \\ 1& x_{n,1} & \cdots & x_{n,p} \end{pmatrix}
\end{equation}$$

 

이때 $L(\beta)$를 $\beta$에 대해서 미분한 후 0으로 놓고 풀면 가중 최소 제곱 회귀 추정량 $\hat{\beta}^{\text{WLS}}$는 다음과 같이 계산할 수 있습니다.

 

$$\hat{\beta}^{\text{WLS}}=(X^tWX)^{-1}X^tWy$$

- 언제 사용하는가? -

가중 최소 제곱법은 다음과 같은 상황에 사용할 수 있습니다.

 

(1) 오차의 등분산성 가정이 의심스러울 때

(2) 이상치(Outlier)에 영향을 덜 받는(Robust) 회귀 모형을 만들고 싶을 때

 

이번 포스팅에서는 상황 (1)에서 가중 최소 제곱법을 적합하는 방법에 대해서 알아보겠습니다. 여기서 오차의 등분산성 가정이 의심스럽다는 것은 오차의 분산이 설명변수에 영향을 받아 관측치에 따라 분산이 달라지는 경우를 말합니다. 아래 그림은 이와 같은 상황을 나타내는 잔차의 형태입니다.

 

 

위와 같은 경우에 가중 최소 제곱 회귀모형을 사용합니다. 

- 가중치를 어떻게 잡아야 하는가? -

만약 분산을 아는 경우(그럴 일은 거의 없지만)에는 즉, 모든 $i$에 대해서 $\sigma_i^2$을 아는 경우에는 $w_i=1/\sigma_i^2$로 설정하여 가중 최소 회귀 추정량을 구합니다.

 

분산을 모르는 경우에는 가중치를 다음과 같이 잡아줍니다.

 

1) 먼저 가중치 없는 최소제곱법을 이용하여 회귀 모형을 구하고 이 모형에 해당하는 잔차 $e_i$를 구합니다.

 

2) 잔차의 절대값 $|e_i|$를 반응 변수, 오차 분산에 영향을 주는 변수 $\tilde{x}_i$를 설명 변수로 하는 가중치 없는 최소 제곱법을 이용하여 회귀 모형을 구합니다. 

 

3) 2)에서 구한 회귀 계수를 $\hat {\eta}_0, \ldots ,\hat {\eta}_p$ 라고 하면 $x_i$에 대응하는 $i$번째 적합 값 $s_i=\hat {\eta}_0+\hat {\eta}_1x_{i1}+\cdots+\hat {\eta}_px_{ip}$를 구하고 $w_i = 1/s_i^2$으로 설정한 후 가중 최소 제곱 법을 이용하여 회귀 모형을 구합니다. 

 

4) 3)에서 구한 회귀 계수와 1)에서 구한 회귀 계수의 차이가 크지 않다면 3)에서 구한 회귀 모형을 최종 모형으로 하고 그렇지 않다면 3)에서 구한 모형을 이용하여 잔차를 계산한 뒤 2) 단계로 넘어갑니다.

- 회귀 계수의 95% 근사 신뢰구간 -

$\hat{\beta}^{\text{WLS}}$의 $j$번째 계수의 95% 근사 신뢰구간은 다음과 같습니다.

$$(\hat{\beta}^{\text{WLS}}_j-s_wt_{0.975, n-p},\;\; \hat{\beta}^{\text{WLS}}_j+s_wt_{0.975, n-p})$$

여기서 $t_{0.975, n-p}$는 자유도가 $n-p$ 인 $t$ 분포의 오른쪽 꼬리 확률이 0.25가 되는 값이고

$$s_w^2=\frac{\sum_{i=1}^nw_ie_i^2}{n-p}$$

입니다.

 

이렇게 구한 신뢰구간은 오차의 분산을 모르는 경우 근사적인 신뢰구간이라고 알려져 있습니다(이런 이유로 95% 근사 신뢰구간이라고 했습니다). 따라서 샘플 수에 따라서 신뢰구간의 정확성이 달라질 수 있지요. 

- 가중 최소 제곱 회귀 모형의 장단점 -

1) 장점

오차의 등분산성을 만족하지 않는 상황에서 가중 최소 제곱 회귀 모형은 어떠한 장점이 있을까요? 설명의 편의를 위해 오차의 분산은 안다고 가정해볼게요. 먼저 반응 변수 $y$와 설명 변수 행렬 $X$의 관계는 다음과 같습니다.

$$y = X\beta + \epsilon\tag{2}$$

여기서 $\epsilon = (\epsilon_1, \ldots, \epsilon_n)^t$이고

$$E(\epsilon)=0 $$ $$ \begin{equation} 
Var(\epsilon)= 
\begin{pmatrix} \sigma_1^2 & 0 & \cdots & 0 \\ 0& \sigma_2^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0& 0 & \cdots & \sigma_n^2 \end{pmatrix} 
\end{equation}. $$

 

그리고 다음의 $n \times n$ 행렬 $W^{1/2}$을 정의합니다.

 

$$\begin{equation} 
W^{1/2}= 
\begin{pmatrix} \sqrt{w_1} & 0 & \cdots & 0 \\ 0& \sqrt{w_2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0& 0 & \cdots & \sqrt{w_n} \end{pmatrix} 
\end{equation}$$

 

여기서 $w_i = 1/\sigma_i^2$입니다. 이제 $W^{1/2}$를 식 (2) 양변 왼쪽에 곱해줍니다.

 

$$W^{1/2}y = W^{1/2}X\beta + W^{1/2}\epsilon \tag{3}$$

 

$y_w=W^{1/2}y$, $ X_w=W^{1/2}X$, $\epsilon_w=W^{1/2}\epsilon$라 하면 식 (3)은 다음과 같이 쓸 수 있습니다.

 

$$y_w = X_w\beta+\epsilon_w\tag{4}$$

 

식 (4)에서 변환된 오차 $\epsilon_w$의 기대값과 분산은 다음과 같습니다.

 

$$E(\epsilon_w) = W^{1/2}E(\epsilon) = 0 \\ Var(\epsilon_w) = W^{1/2}Var(\epsilon)W^{1/2} = I$$

 

즉, 식 (4)는 오차의 기대값이 0이며 $Var(\epsilon_{wi})=1$이므로, 등분산성을 만족하는 것을 알 수 있습니다. 식 (4)에서 최소 제곱법을 이용하여 구한 회귀 추정량 $\hat{\beta}_w$는 다음과 같습니다.

 

$$\hat{\beta}_w = (X_w^tX_w)^{-1}X_w^ty_w = (X^tWX)^{-1}X^tWy$$

 

이 결과를 통해 가중 최소 제곱 회귀 추정량은 식 (4)에서와 같이 오차의 등분산성을 만족하는 변환된 데이터의 일반 최소 제곱 회귀 추정량과 같아집니다. 또한 가중 최소 제곱 추정량 $\hat{\beta}_w$은 불편추정량이며 따라서 $\hat{\beta}_w$은 Gauss-Markov 정리에 의하여 불편 선형 추정량 중에서 최소 분산을 갖는 최적(optimal) 추정량이 됩니다. 이러한 장점이 있어서 오차의 등분산성이 만족되지 않는 경우 가중 최소 제곱 회귀 모형을 사용하는 것입니다.

 

2) 단점

오차의 등분산성 가정이 위배된 경우 가중 최소 제곱 추정량은 일반적인 최소 제곱 추정량과 같이 이상치(Outlier)에 민감하다는 단점이 있습니다. 오차의 분산을 모르는 경우 이를 잔차의 절대값과 설명변수로 이루어진 회귀 모형으로 추정을 하는데요. 데이터의 개수가 적은 경우 추정의 정확성이 떨어져 가중 최소 제곱 회귀 추정량을 신뢰할 수 없게 됩니다. 또한 가중 최소 제곱 회귀 추정량의 신뢰구간은 근사적으로 구하는 것이기 때문에 데이터의 개수가 적을 경우 신뢰구간의 정확성도 떨어지게 됩니다.

반응형

---

   2. 가중 최소 제곱법(Weighted Least Square) 적합하기 with Python

먼저 데이터를 다운받아주세요.

 

blood_pressure.csv

0.00MB

blood_pressure_description.txt

0.00MB

 

다음으로 필요한 모듈을 임포트하고 데이터를 불러옵니다.

 

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

from scipy.stats import t
from statsmodels.formula.api import ols
from statsmodels.formula.api import wls

df = pd.read_csv('./blood_pressure.csv') ## 데이터 불러오기

 

 

먼저 오차의 등분산성을 확인하기 위해 잔차도를 그려볼 거예요. 아래 코드는 설명 변수(Age)에 따른 잔차의 산포도를 그린 것입니다.

 

fit = ols('DBP ~ Age',data=df).fit() ## 단순선형회귀모형 적합

fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
ax = fig.add_subplot()
ax.scatter(df['Age'],fit.resid) ## 잔차도 출력
ax.axhline(y=0, color='red')
plt.xlabel('Age', fontsize=15)
plt.ylabel('Residual', fontsize=15)
plt.show()

 

 

위의 잔차도를 보시면 오차의 등분산성이 만족되지 않음을 알 수 있습니다. 가중 최소 제곱법을 이용하여 회귀 모형을 구하기 전에 (가중치 없는) 최소 제곱 회귀 추정량을 살펴봅시다.

 

print(fit.params.Age)
print(fit.params.Intercept)

 

 

따라서 최소 제곱 회귀 모형은 다음과 같다는 것을 알 수 있습니다.

$$\widehat{\text{DBP}}^{\text{LS}} = 56.1569 + 0.5800 \text{AGE}\tag{5}$$

 

이제 가중 최소 제곱 회귀 모형을 구해봅시다.

 

df['Abs_Residual'] = abs(fit.resid) ## 잔차 절대값
std_fit = ols('Abs_Residual~Age',data=df).fit() ## 잔차 절대값과 연령을 이용하여 회귀 모형 적합

 

line 1~2

앞서 구한 회귀 모형을 이용한 잔차(fit.resid)를 절대값으로 변환합니다(line 1). 이 값을 반응 변수, Age를 설명변수로 하는 회귀 모형을 적합합니다(line 2).

 

std_fit.summary()

 

line 1

회귀 계수가 유의미한지 확인해보기 위하여 summary를 이용했습니다.

 

 

코드를 실행하면 다음과 같이 모형에 대한 정보가 출력됩니다. 빨간색 네모 박스를 보시면 Age에 대응하는 회귀 계수의 $p$-value가 굉장히 작다는 것을 확인할 수 있고 따라서 회귀 계수가 유의미하다는 것을 알 수 있습니다. 

 

이제 가중 최소 제곱 회귀 모형을 적합해보겠습니다. 아래 코드는 가중 최소 제곱 회귀 모형을 구현한 것입니다.

 

p = 2 ## intercept, slope
y = df['DBP'] ## Response
weight = 1/np.square(std_fit.fittedvalues)
W = np.diag(weight) ## Weight Diagonal Matrix

## Model Matrix
X = pd.DataFrame()
X['Intercept'] = [1]*len(df)
X['Age'] = df['Age']
X = np.array(X)

XtWX = np.matmul(np.matmul(X.transpose(),W),X)
XtWX_inv = np.linalg.inv(XtWX)
XtWy = np.matmul(np.matmul(X.transpose(),W),y)
weighted_coef = np.matmul(XtWX_inv,XtWy) ## weighted least square estimator of slope

 

line 3~4

$s_i$값을 구하고 제곱의 역수를 가중치로 잡았습니다(line 3). 그리고 이 가중치를 대각 원소로 하는 대각 행렬 $W$를 만들었습니다.

 

line 7~15

가중 최소 제곱 회귀 추정량 $\hat{\beta}^{\text{WLS}}=(X^tWX)^{-1}X^tWy$를 계산합니다.

 

$\hat{\beta}^{\text{WLS}}$을 확인해보겠습니다.

 

 

위 그림을 보면 가중 최소 제곱 회귀 모형은 다음과 같다는 것을 알 수 있습니다.

$$\widehat{\text{DBP}}^{\text{WLS}} = 55.5658 + 0.5963 \text{AGE}\tag{6}$$

 

가중 최소 제곱 회귀 모형 (6)과 최소 제곱 회귀 모형 (5)을 비교해보았을 때 회귀 계수는 그렇게 차이가 나지 않는 것으로 판단됩니다. 따라서 (6)을 최종 모형으로 선택합니다.

 

다음으로 $\hat{\beta}^{\text{WLS}}$의 95% 근사 신뢰구간을 구해보았습니다.

 

## confidence interval of slope
MSE_w = np.sum(weight*np.square(fit.resid))/(len(df)-p)
coef_cov_matrix = MSE_w*XtWX_inv ## estimated covariance matrix of weighted least square coefficient

se_slope = np.sqrt(coef_cov_matrix[1][1]) ## estimated standard deviantion of weighted least square coefficient 
critical_value = t.ppf(0.975, len(df)-p)
lower_limit = weighted_coef[1] - critical_value*se_slope 
upper_limit = weighted_coef[1] + critical_value*se_slope

 

 

신뢰구간이 0을 포함하지 않으므로 가중 최소 제곱 회귀 추정량 $\hat{\beta}^{\text{WLS}}$은 유의미하다고 볼 수 있습니다.

 

이제 $\hat{\beta}^{\text{WLS}}$을 이용하여 다시 잔차도를 그려보겠습니다.

 

## 가중치를 이용하여 변수 변환
sqrt_weight = np.sqrt(weight)
sqrt_W = np.diag(sqrt_weight)

y_w = np.matmul(sqrt_W,y)
X_w = np.matmul(sqrt_W,X)
weighted_residual = y_w - np.matmul(X_w,weighted_coef)
weighted_AGE = X_w[:,1]

## 잔차도 출력
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
ax = fig.add_subplot()
ax.scatter(weighted_AGE,weighted_residual) 
ax.axhline(y=0, color='red')
plt.xlabel('weighted_AGE', fontsize=15)
plt.ylabel('weighted_residual', fontsize=15)
plt.show()

 

line 2~8

$\hat{\beta}^{\text{WLS}}$을 이용하여 잔차도를 그리기 위해서는 식 (3)과 같이 가중치를 이용한 변수 변환을 해주어야 합니다.

 

위 코드를 실행해보세요.

 

 

보시는 바와 같이 설명변수에 대해서 잔차의 특정 패턴이 없어졌음을 알 수 있고 이에 따라 오차의 등분산성이 만족된다고 판단할 수 있습니다.

반응형

---

   3. statsmodels 모듈을 이용하여 가중 최소 제곱 모형 적합하기

이번에는 statsmodels 모듈을 이용하여 가중 최소 제곱 회귀 모형을 적용하는 방법에 대해서 알아보겠습니다.

 

df = pd.read_csv('./blood_pressure.csv') ## 데이터 불러오기
fit = ols('DBP ~ Age',data=df).fit() ## 단순선형회귀모형 적합

df['Abs_Residual'] = abs(fit.resid)
std_fit = ols('Abs_Residual~Age',data=df).fit() ## 잔차의 절대값 ~ 연령
wls_fit = wls('DBP ~ Age',data=df, weights = 1/np.square(std_fit.fittedvalues)).fit()

 

line 6

statsmodels에서 제공하는 wls를 이용하여 가중 최소 제곱 회귀 모형을 적합합니다. 모형 표현식은 ols와 동일하며 가중치를 weights인자에 넣어주기만 하면 됩니다.

 

summary를 이용하여 결과를 확인해보겠습니다.

 

 

회귀 계수 추정값과 신뢰구간을 보시면 위에서 구한 값과 동일하다는 것을 알 수 있습니다.

---

이번 포스팅에서는 가중 최소 제곱법이 무엇인지 적합하는 방법은 어떻게 되는지 알아보았습니다. 모형 적합의 경우 그냥 간단하게 statsmodels를 이용하면 됩니다. 하지만 모형 적합 과정을 알아야 결과를 해석하고 이해하는데 도움이 되므로 적합 알고리즘을 이해하셨다면 직접 구현도 해보시는 것을 추천합니다.

 

지금까지 꽁냥이의 글 읽어주셔서 감사합니다. 안녕히 계세요.

---