Multiple Linear Regression(다중 선형 회귀)


  • 본문의 일부는 Data Science from scratch(밑바닥부터 시작하는 데이터 과학)을 참고했습니다.
  • 본문의 소스코드의 일부는 Joel Grus의 Github에서 Unlicensed 라이선스로 배포되고 있습니다.

Features, 특성

Supervised Learning(교사학습) 모델에서 학습을 통해 찾고자 하는 목표값(y)을 위해 현실에서는 엄청나게 많은 x를 이용할 수 있다. 데이터 과학에서 x는 feature(혹은 attribute, variable, predictor 등으로도 불린다. 한국어로는 특성, 속성 등으로 번역할 수 있다.)라고 불리며, 실무에서 데이터를 다룰 경우 처럼 feature가 많은 경우 모든 feature가 모델의 성능을 끌어올리는 데에 도움이 될 수도 있지만 그렇지 않을 수도 있다.

The initial set of raw features can be redundant and too large to be managed. Therefore, a preliminary step in many applications of machine learning and pattern recognition consists of selecting a subset of features, or constructing a new and reduced set of features to facilitate learning, and to improve generalization and interpretability. - Wikipedia, Feature selection

Feature들의 초기집합은 관리하기에 너무 크고 많을 수 있다. 그렇기 때문에 많은 기계학습, 패턴 인식의 예비 단계에서는 feature를 선택하거나 양을 줄여 새로운 feature를 만들어 학습과 일반화, 해석을 용이하게 만들 수 있도록 한다.

모델에 적합한 feature를 찾기 위한(Feature selection) 여러가지 방법론이 존재하며, 그 중 상관관계 특성선택에서는 ‘Good feature subsets contain features highly correlated with the classification, yet uncorrelated to each other(좋은 feature 서브셋은 분류와 상관관계가 아주 높고 서로간의 상관관계가 적다.)’라는 가설을 갖기도 한다. 또 feature를 통해 학습하는 것 뿐만이 아니라 feature 자체를 학습하는 feature learning(특성학습)도 존재한다.

다중 회귀 분석(Multiple linear regression)에서는 단순 회귀 분석에서보다 많은 feature를 통해 모델의 성능을 끌어올리는 것을 목표로 한다.

동기부여

사용자들의 하루 평균 근무 시간, 박사 학위 취득 여부, 친구의 수라는 특성을 이용해 사이트 체류시간 예측의 성능을 높여보고자 한다. 먼저 독립 변수를 사용하는 선형 모델을 시험해보자.

시간(분) = \(\alpha + \beta_1 * (친구 수) + \beta_2 * (근무 시간) + \beta_3 * (박사 학위 취득 여부) + \epsilon\)

박사 학위 취득 여부는 숫자 데이터가 아니지만 참, 거짓(0과 1) 여부로 가변수(dummy variable)를 만들어 표현할 수 있다.

모델

단순 회귀 분석에서는 아래와 같은 모델을 다뤘다.

\[y_i = \alpha + \beta x_i + \epsilon_i\]

각 입력값 \(x_i\)을 숫자 하나가 아닌 k개의 숫자인 \(x_i,...,x_k\)라고 하면 다중 회귀(multiple regression)모델은 다음과 같아진다.

\[y_i = \alpha + \beta_i x_{i1} +...+ \beta_k x_k + \epsilon-i\]

다중 회귀 분석에서 파라미터 벡터를 보통 \(\beta\)라고 부른다. 상수항 \(\alpha\)까지 덧붙이려면 각 데이터 x_i의 앞부분에 1을 덧붙이면 된다.

beta = [alpha, beta_1, ..., beta_k]

그러면 각 데이터는 다음과 같다.

x_i = [1, x_i1, ..., x_ik]

이렇게 하면 모델을 다음처럼 나타낼 수 있다.

def predict(x_i, beta):
  """각 x_i의 첫 항목은 1이라고 가정"""
  return dot(x_i, beta)

""" 독립 변수 x는 각각 다음과 같은 벡터의 열로 표현할 수 있다. """

[1,   # 상수항
49,   # 친구의 수
4,    # 하루 근무 시간
0]    # 박사 학위 취득 여부

최소자승법에 대한 추가 가정

이 모델이 의미를 갖기 위해 몇 가지 추가적 가정이 필요하다.
먼저 x의 열은 서로 일차독립해야 한다. 일차독립이란 어떤 벡터도 다른 벡터의 선형결합으로 만들어질 수 없다는 걸 의미한다. 이 가정이 성립하지 않는다면 베타를 추정할 수 없다.
극단적인 예로 num_friends(친구 수)와 같은 num_acquaintances(지인 수)가 데이터에 추가된다고 가정하면 어떠한 beta를 사용해도 num_friends 계수에 임의의 값을 더하고 num_acquaintances 계수에서 같은 값을 빼주면 모델의 예측값은 변하지 않는다. 즉 num_friends의 계수를 계산할 수 없다는 뜻이다. (보통 이 가정이 위배됐는지 확인하기란 쉽지 않다.)
두 번째로 중요한 가정은 x의 모든 열은 오류 \(\epsilon\)과 아무런 상관 관계가 없다는 것이다. 만약 이에 위배되는 경우 아예 잘못된 beta가 추정된다. 예컨데 단순 회귀 모델을 살펴보면 친구 수가 한 명 증가할 때 사용자의 하루 평균 사이트 체류 시간이 0.9분 증가한다고 예측됐다.
다음과 같은 경우가 있을 수도 있다.

  • 근무 시간이 더 긴 사람은 더 적은 시간을 사이트에서 보낼 것이다.
  • 친구 수가 많은 사람일수록 근무 시간이 더 길다.

그렇다면 실제 모델은 다음과 같고 근무 시간과 친구 수는 양의 상관관계를 지닌다고 가정하자.

사이트에서 보내는 시간(분) = \(\alpha + \beta_1 친구 수 + \beta_2 근무 시간 + \epsilon\)

만약 다음과 같은 단순 회귀 분석 모델에서 오류를 최소화시키면 \(\beta_i\)는 과소평가될 것이다.

사이트에서 보내는 시간(분) = \(\alpha+\beta_1 친구 수 + \epsilon\)

만약 실제 모델의 \(\beta_1\)을 사용하는 단순 회귀 모델로 예측 성능을 평가해 보자. \(\beta_2\) > 0이지만 모델에 포함시키지 않았기에 근무 시간이 긴 사용자의 예측값은 너무 작게 계산될 것이고 근무 시간이 큰 사용자의 예측값은 너무 크게 계산될 것이다. 또 근무 시간과 친구의 수는 양의 상관관계를 지녀 친구 수가 많은 사용자의 예측값은 너무 작게 계산될 것이고 친구 수가 적은 사용자의 예측값은 너무 크게 계산될 것이다.
단순 회귀 모델의 오류를 줄이기 위해선 추정된 \(\beta_1\)을 줄여야 한다. 즉 오류를 최소화하는 \(\beta_1\)은 실제 값보다 작아진다는 뜻이다. 결국 단순 회귀 모델은 \(\beta_1\)을 과소평가하도록 편향된다. 보통 이렇게 독립 변수와 오류 사이에 상관관계가 존재할 경우 최소자승법으로 만들어지는 모델은 편향된 \(\beta\)를 추정해준다.

모델 학습하기

단순 회귀 모델처럼 sum(오류**2)값을 최소화 해주는 beta를 찾아야 한다. 이 때 경사하강법을 사용한다. 먼저 최소화할 오류 함수를 만들자. 일단 SGD(Stochastic Gradient Descent)를 사용하기 위해 예측값 하나에 대한 오류의 제곱 값이 필요하다.

import random
from linear_algebra import dot
from simple_linear_regression import total_sum_of_squares
from gradient_descent import minimize_stochastic

def predict(x_i, beta):
    """각 x_i의 첫 항목은 1이라고 가정"""
    return dot(x_i, beta)

def error(x_i, y_i, beta):
    return y_i - predict(x_i, beta)

def squaredErrors(x_i, y_i, beta):
    return error(x_i, y_i, beta) ** 2

def squaredErrorGradient(x_i, y_i, beta):
    """i번째 오류 제곱 값의 beta에 대한 기울기"""
    return [-2 * x_ij * error(x_i, y_i, beta) for x_ij in x_i]

def estimateBeta(x, y):
    betaInitial = [random.random() for x_i in x[0]]
    return minimize_stochastic(squaredErrors,
                               squaredErrorGradient,
                               x, y,
                               betaInitial,
                               0.001)
x = [[some_list_with_four_arguments_filled_in_each]]
dailyMinutes = [list_with_numbers]

beta = estimateBeta(x, y)
>> # 아래와 같은 모델이 만들어진다.
[30.619881701311712,
 0.9702056472470465,
 -1.8671913880379478,
 0.9163711597955347]

 """ 분 = 30.619 + 0.97 친구 수 - 1.86 근무 시간 + 0.916 박사 학위 여부"""

모델 해석하기

모델의 계수는 다른 모든 조건이 동일할 때 해당 항목의 영향력을 뜻한다. 다른 모든 게 동일할 때 친구 수가 한 명 증가하면 사용자가 하루 평균 사이트를 이용하는 시간이 1분 늘어난다. 또 다른 조건이 동일할 때 근무 시간이 한 시간 늘어나면 사이트 이용 시간이 2분 줄어든다. 박사 학위가 있는 사람은 다른 조건이 동일할 때 사이트를 1분 더 이용한다.

이런 해석으론 변수 간의 관계를 직접 설명할 수 없다. 예컨데 근무 시간이 서로 다른 사용자들의 친구 수가 다를 수도 있다. 이 모델은 이런 관계를 잡아내지 못한다. 이런 문제를 해결하기 위해 친구 수와 근무 시간을 곱해 새로운 변수를 만들면 새 변수를 통해 친구 수의 증가에 따른 근무 시간의 계수를 증가/감소 시킬 수 있다.

아니면 친구 수가 증가할수록 사이트에서 보내는 시간이 일정 수준까지 증가하고 그 후엔 감소할 수도 있다. (친구가 너무 많은 게 부담스러웠던 걸까?) 이런 현상을 잡기 위해 친구 수를 제곱한 값을 변수로 사용할 수 있다. 변수가 점점 추가되면 각 계수가 유의미한지 확인해야 한다. 변수끼리의 곱, 변수의 log값, 변수의 제곱 값 등 수 많은 변수를 추가할 수 있기 때문이다. 중요한 건 어떤 가설을 생각해내고 그에 맞는 변수들을 테스트해보는 것이다.

적합성(Goodness of fit)

모델의 R 제곱 값(에러의 제곱의 합 / 총 제곱의 합)을 측정해보자. 이전에 사용한 모델보다 발전된 모습을 볼 수 있다.

def multipleRSquared(x, y, beta):
  sumOfSquaredErrors = sum(error(x_i, y_i, beta) ** 2
                          for x_i, y_i in zip(x,y))
  return 1.0 - sumOfSquaredErrors / total_sum_of_squares(y)

multipleRSquared(x, daily_minutes, beta)
>> 0.6800074955952597

하지만 사실 회귀 분석 모델에 새 변수를 추가하면 R 제곱 값이 어쩔 수 없이 증가할 수밖에 없다. 단순 회귀 분석 모델은 근무 시간과 박사 학위 취득 여부가 0인 다중 회귀 모델과 동일하다. 즉 최적의 다중 회귀 분석 모델은 언제나 단순 회귀 모델보다 작은 오류를 갖는다.

그렇기 때문에 다중 회귀 분석 모델에서는 각 계수의 표준 오차를 살펴봐야 한다. 표준 오차는 추정된 \(\beta_1\)의 계수가 얼마나 확실한지 보여준다. 모델은 주어진 데이터에 적합(fit)할 수도 있으나, 몇몇 독립 변수 간에 어떤 상관관계가 있다면 이 변수들의 계수는 무의미할 것이다.
오차를 측정하기 위해선 각 오류 \(\epsilon_1\)은 독립이며 평균은 0이고 표준편차는 \(\sigma\)인 정규분포의 확률변수라는 가정이 필요하다. 선형대수를 통해 각 계수들의 표준 오차를 계산할 수 있는데, 표준 오차가 클수록 해당 계수는 무의미해진다.

bootstrap

(여기서 말할 bootstrap은 트위터에서 재창한 어떤 웹 UI 프레임워크와는 별 관계가 없다)
만약에 알 수 없는 분포에서 생성된 표본 데이터가 주어졌다고 상상해보자.

data = getSample(num_points=n)

중앙값은 분포 자체의 중앙값에 대한 추정치로 사용할 수 있다. 하지만 만약 표본 데이터가 모두 100 근처에 위치한다면 중앙값 또한 100 근처에 위치할 것이다. 만약 데이터의 반이 0 근처에 위치하고, 나머지 절반이 100 근처에 위치한다면 추정된 중앙값은 신뢰하기 힘들다.
만약 새로운 표본을 계속 얻을 수 있다면, 각 표본의 중앙값을 계산해 보며 값들의 분포를 살펴볼 수 있다. 하지만 보통 새로운 표본을 계속 받을 수 없기 때문에, bootstrap이란 새로운 표본 데이터를 만드는 기법을 사용할 수 있다.

bootstrap은 기존의 데이터에서 중복 가능한 재추출을 통해 새 데이터의 항목을 생성한다. 그리고 그 인공 데이터로 중앙값을 계산해 볼 수 있다. 데이터를 만들어서 실험해보자.

def bootstrapSample(data):
    """무작위 값을 len(data)만큼 재추출"""
    return [random.choice(data) for _ in data]

def bootstrapStatistics(data, statsFn, numSamples):
    return [statsFn(bootstrapSample(data)) for _ in range(numSamples)]

closeTo100 = ( [ 99.5 + random.random() for _ in range(101)] )
farFrom100 = ( [99.5 + random.random()] +
               [random.random() for _ in range(50)] +
               [200 + random.random() for _ in range(50)] )

print("Close to 100 : ",bootstrapStatistics(closeTo100, median, 10))
print("")
print("Far from 100 : ",bootstrapStatistics(farFrom100, median, 10))
>>
Close to 100 :  [99.8514476396684, 100.01804048374089, 99.96781493910923, 99.98715682969906, 100.0208555789534, 99.9962802757804, 99.98715682969906, 100.00779191205707, 100.01804048374089, 99.90836334359076]

Far from 100 :  [0.8854546218754292, 200.05648615607967, 200.03521556425756, 200.05648615607967, 0.9779031140920793, 0.9008187232236743, 200.05648615607967, 200.0370302969729, 0.9303349466491837, 0.971366967742836]

standard_deviation(bootstrapStatistics(farFrom100, median, 10))
>>
99.00177259865933
standard_deviation(bootstrapStatistics(closeTo100, median, 10))
>>
0.05488572073464672

계수의 표준 오차

계수의 표준 오차를 추정할 때도 bootstrap을 활용할 수 있다. 주어진 데이터를 bootstrap 함수에 넣어 다수의 bootstrap 데이터를 생성해 각각의 데이터에서 beta를 추정할 수 있다. 만약 모든 boostrap 데이터에서 친구 수에 해당하는 독립 변수의 계수가 그다지 변하지 않는다면 추정된 계수는 신뢰할 수 있다. 반대로 계수가 많이 변한다면 추정된 계수는 신뢰할 수 없다.
유의할 점은 x와 y를 zip으로 묶어야 하는 점이다.

def estimateSampleBeta(sample):
    """sample은 (x_i,y_i)의 list"""
    x_sample, y_sample = zip(*sample)    # zip을 풀어주는 마법
    return estimateBeta(x_sample, y_sample)

bootstrapBetas = bootstrapStatistics(list(zip(x, daily_minutes)),
                                    estimateSampleBeta,
                                    100)    # 시간이 좀 걸린다.

bootstrapStandardErrors = [ standard_deviation([beta[i] for beta in bootstrapBetas]) for i in range(4)]

bootstrapStandardErrors
>>

[0.953551702104508,   # 상수,        실제 오류 = 1.19
 0.06288763616183773, # 친구의 수,   실제 오류 = 0.080
 0.11722269488203318, # 근무 시간,   실제 오류 = 0.127
 0.8591786495949066]  # 박사 학위,   실제 오류 = 0.998  

이제 \(\beta_1\)는 0일까 같은 가설을 검증할 수 있다. \(\beta_1 = 0\)이란 귀무가설을 확인하기 전에 분포에 대한 기본적인 가정과 함께 다음과 같은 통계치를 계산할 수 있다.

\[t_j=\hat{\beta_j}/\hat{\sigma_j}\]

(\(\hat{\beta_j}\)는 샘플링을 통해 계산한 \(\beta\)를 의미한다.)
\(\beta_j\)의 추정치를 표준 오차의 추정치로 나눈 값은 n-k 자유도를 지닌 t 분포를 따른다.
만약 t 분포의 누적분포함수를 계산해 주는 students_t_cdf 함수가 있었다면 각 계수의 p-value를 계산해 가설을 검증해 볼 수 있었을 것이다. 하지만 아쉽게도 그런 함수를 만든 적은 없다.
그러나 자유도가 높아질수록 t 분포는 표준정규분포와 비슷해진다. n이 k보다 훨씬 큰 경우 normal_cdf를 사용할 수 있다.

def pValue(betaHatJ, sigmaHatJ):
    if betaHatJ > 0:
        return 2 * (1 - normal_cdf(betaHatJ / sigmaHatJ))
    else:
        return 2 * normal_cdf(betaHatJ / sigmaHatJ)

print(beta,"\n",bootstrapStandardErrors)
>>
[30.619881701311712, 0.9702056472470465, -1.8671913880379478, 0.9163711597955347]
[0.953551702104508, 0.06288763616183773, 0.11722269488203318, 0.8591786495949066]

print(pValue(30.619, 0.970)) # ~0
print(pValue(0.97, 0.062))   # ~0
print(pValue(-1.867, 0.117)) # ~0
print(pValue(0.916, 0.859))  # 0.28626264064150586

(나중에는 정확한 t분포의 값과 표준 오차 계산을 위해 다른 라이브러리를 사용하도록 하자.) 대부분의 계수에서 p-value가 굉장히 작기 때문에 이 계수들은 0이 아닌 것으로 검증됐다. 하지만 박사 학위 취득 여부의 p-value가 0보다 큰 값으로 계산되어 이 계수의 의미가 없을 수도 있다는 것이 암시된다.

Regularization

실제 데이터를 분석할 때는 변수가 굉장히 많은 데이터에 회귀 모델을 적용하는 경우가 많다. 그러나 변수가 너무 많으면 다양한 문제가 발생할 수 있다. 첫 번째로 변수가 많아지면 모델은 당연히 학습 데이터에 오버피팅할 것이다. 또 0이 아닌 계수가 많아질수록 모델을 해석하는 게 어려워진다. 만약 어떤 현상을 설명하는 게 목표라면 수백 개의 변수로 모델을 만드는 것보다 세 개 정도의 변수로 작은 모델을 만드는 게 나을 것이다.
Regularization(정규화, Normalization은 데이터 범위를 유사하게 맞춰주는 일을 말한다.)은 beta가 커질수록 해당 모델에게 패널티를 주는 방법이다. 오류와 패널티를 동시에 최소화하면 최적의 모델을 만들 수 있다. 패널티를 강조할수록 값이 큰 계수에 더 큰 제한이 걸린다.

ridge regression의 경우 beta_i를 제곱한 값의 합에 비례하는 패널티를 추가한다. 하지만 상수에는 패널티를 주지 않는다.

# alpha는 패널티의 강도를 조절하는 하이퍼파라미터이다.
# 보통은 lambda를 쓰지만 이는 파이썬에서 이미 사용중이다.

def ridge_penalty(beta, alpha):
    return alpha * dot(beta[1:], beta[1:])

def squared_error_ridge(x_i, y_i, beta, alpha):
    """beta를 사용할 때 오류와 패널티의 합 추정"""
    return error(x_i,y_i, beta) ** 2 + ridge_penalty(beta,alpha)

def ridge_penalty_gradient(beta, alpha):
    return [0] + [2 * alpha * beta_j for beta_j in beta[1:]]

# 경사하강법을 통해 학습시킬 수 있다.
def squared_error_ridge_gradient(x_i, y_i, beta, alpha):
    return vector_add(squaredErrorGradient(x_i, y_i, beta),
                     ridge_penalty_gradient(beta, alpha))

def estimate_beta_ridge(x,y,alpha):
    """패널티가 alpha인 리지 회귀를 경사하강법으로 학습"""
    beta_initial = [random.random() for x_i in x[0]]
    return minimize_stochastic(partial(squared_error_ridge, alpha=alpha),
                               partial(squared_error_ridge_gradient,
                                      alpha=alpha),
                                      x,y,
                                      beta_initial,
                                      0.001)
                                      random.seed(0)

beta_0 = estimate_beta_ridge(x, daily_minutes, alpha=0.0)
print(multipleRSquared(x, daily_minutes, beta_0))
print(beta_0)
beta_0 = estimate_beta_ridge(x, daily_minutes, alpha=0.01)
print(multipleRSquared(x, daily_minutes, beta_0))
print(beta_0)
beta_0 = estimate_beta_ridge(x, daily_minutes, alpha=0.1)
print(multipleRSquared(x, daily_minutes, beta_0))
print(beta_0)
beta_0 = estimate_beta_ridge(x, daily_minutes, alpha=10)
print(multipleRSquared(x, daily_minutes, beta_0))
print(beta_0)
>>
0.6800074955952597
[30.619881701311712, 0.9702056472470465, -1.8671913880379478, 0.9163711597955347]
0.680010213297079
[30.55985204967343, 0.9730655363505671, -1.8624424625144256, 0.9317665551046306]
0.6797276241305292
[30.894860179735474, 0.9490275238632391, -1.8501720889216575, 0.5325129720515789]
0.5745539861941786
[28.43248665430619, 0.7246794348835753, -0.9225528545940292, -0.01951794239121216]    

패널티가 증가하면 증가할수록 박사 학위 여부 변수가 사라져가는 걸 확인할 수 있다. 이는 박사 학위 여부라는 feature 자체가 그다지 유효하지 않다는 걸 의미한다.

다른 형태의 패널티를 활용하는 lasso regression도 있다.

def lasso_panelty(beta, alpha):
  return alpha * sum(abs(beta_i) for beta_i in beta[1:])

리지 회귀의 패널티는 총 계수의 합을 줄여 주지만 라쏘 회귀의 패널티는 모든 계수를 최대한 0으로 만들어 주며 더 희소한(sparse) 모델을 학습하게 해준다. 하지만 경사하강법을 통해 학습을 시킬 수 없다.