5.8 비선형 회귀

지금까지 이 장에서 선형 관계를 가정하는 것이 종종 적절하긴 했지만, 비선형 형태가 더 적당한 경우가 많이 있습니다. 이 절에서는 내용을 좀 더 쉽게 설명하기 위해, 단 하나의 예측변수(predictor variable) \(x\)만 있다고 가정합니다.

비선형적인 관계를 모델링하는 가장 단순한 방법은 회귀 모델을 추정하기 전에 목표 예상변수(forecast variable) \(y\)와 예측변수(predictor variable) \(x\)를 변환하는 것입니다(또는 \(y\)나 \(x\) 중의 하나를). 여기에서 비선형적인 함수 형태를 얻긴 하지만, 모델은 매개변수에 대해 여전히 선형적입니다. 가장 흔히 사용하는 변환은 (자연) 로그입니다(3.2 절을 보시길 바랍니다).

로그-로그 함수 형태는 다음과 같습니다. \[ \log y=\beta_0+\beta_1 \log x +\varepsilon. \] 이 모델에서, 기울기 \(\beta_1\)는 탄력성으로 해석할 수 있고, \(\beta_1\)은 \(x\)가 \(1\%\) 증가할 때 \(y\)의 평균 백분율 변화입니다. 다른 유용한 형태들도 적을 수 있습니다. 로그-선형 형태는 목표 예상변수(forecast variable)만 변환하여 적고 선형-로그 형태는 예측변수(predictor variable)를 변환하여 얻습니다.

변수를 로그 변환하기 위해서는, 모든 관측된 값이 0보다 커야한다는 것을 떠올려봅시다. 변수 \(x\)가 0인 값을 가지고 있는 경우에는 \(\log(x+1)\)와 같은 변환을 사용합니다. 즉, 변수 값에 1을 더하여 로그를 취합니다. 이렇게 하여 0 때문에 생기는 문제를 피하고 로그를 취하는 것과 비슷한 효과를 냅니다. 이 경우도 변환된 눈금에서 0이 남아있는 원래의 눈금에서의 0 값 때문에 생기는 부작용을 가지고 있습니다.

단순하게 데이터를 변환하는 것이 적절하지 않을 경우가 있기에 더욱 일반적으로 적는 것이 필요할 수 있습니다. 그러면 사용할 모델은 다음과 같습니다. \[ y=f(x) +\varepsilon \] 여기에서 \(f\)는 비선형 함수입니다. 표준 (선형) 회귀에서, \(f(x)=\beta_{0} + \beta_{1} x\)입니다. 따라나오는 비선형 회귀 형태에서는 \(f\)를 단순히 로그나 다른 변환에 비해 더욱 유연한 \(x\)에 관한 비선형 함수로 둡니다.

가장 단순하게 적는 방법 중의 하나는 \(f\)를 조각별 선형(piecewise linear)으로 만드는 것입니다. 즉, \(f\)의 기울기가 변할 수 있는 점을 도입합니다. 이러한 점을 “매듭(knots)”으로 부릅니다. \(x_{1,t}=x\)로 두고 변수 \(x_{2,t}\)을 도입하여 다음과 같이 적을 수 있습니다. \[\begin{align*} x_{2,t} = (x-c)_+ &= \left\{ \begin{array}{ll} 0 & x < c\\ (x-c) & x \ge c \end{array}\right. \end{align*}\] 기호 \((x-c)_+\)는 양수일 때 \(x-c\) 값을 의미하고 이외에는 0입니다. 이렇게 하여 \(c\) 점에서 기울기가 구부러지도록합니다. 위의 형태에서 추가적인 변수를 더하여 관계식에서 구부러진 부분을 더 넣을 수 있습니다.

이 예제에서는 \(x=t\)로 두고 조각별 선형 추세를 시계열에 맞춥니다.

이런 방식으로 구성한 조각별 선형 관계는 회귀 스플라인(regression splines)의 특수한 경우입니다. 일반적으로, 선형 회귀 스플라인은 아래 관계식을 이용하여 얻습니다. \[ x_{1}= x \quad x_{2} = (x-c_{1})_+ \quad\dots\quad x_{k} = (x-c_{k-1})_+ \] 여기에서 \(c_{1},\dots,c_{k-1}\)은 매듭(선이 휠 수 있는 점)입니다. 매듭(knot)의 수(\(k-1\))를 고르는 것과 어디에 두어야 하는 지는 어려울 수 있고 약간 임의적일 수 있습니다. 자동으로 매듭을 고르는 알고리즘 몇몇을 몇몇 소프트웨어에서 사용할 수 있습니다만, 아직까지 널리 사용되는 것은 아닙니다.

조각별 선 대신 단계적 3차식으로 더 매끄러운 결과를 얻습니다. 이 경우는 제한조건이 걸려있어서 연속적이고(서로 연결되어있는) 매끄럽습니다(따라서 조각별 선형 스플라인에서 나타나는 갑작스런 방향 변화가 없습니다.) 일반적으로 3차(cubic) 회귀 스플라인은 다음과 같이 적습니다. \[ x_{1}= x \quad x_{2}=x^2 \quad x_3=x^3 \quad x_4 = (x-c_{1})_+ \quad\dots\quad x_{k} = (x-c_{k-3})_+. \] 3차 스플라인이 보통 데이터에 더 잘 맞는 결과를 냅니다. 하지만, \(x\)가 과거 데이터의 범위를 넘어갈 때 \(y\) 값 예측값을 신뢰할 수 없게 됩니다.

비선형 추세로 예측하기

5.4 절에서 \(x=t\)로 두고 선형 추세를 시계열에 맞추는 것을 소개했습니다. 비선형 추세를 맞추는 가장 간단한 방법은 다음과 같이 적어서 얻는 2차 이상의 추세를 사용하는 것입니다. \[ x_{1,t} =t,\quad x_{2,t}=t^2,\quad \dots. \] 하지만, 예측 작업에서 2차 이상의 추세를 사용하는 것을 권장하지는 않습니다. 외삽(extrapolate)하는 경우에, 결과 예측값이 종종 비현실적으로 나타납니다.

더 나은 접근 방식은 위에서 소개한 것을 조각별로 적어서 사용하고, 시간에 따라 어떤 점에서 조각별 선형 추세가 휘게 맞추는 것입니다. 이것을 선형 조각으로 구성된 비선형 추세로 생각할 수 있습니다. 시간 \(\tau\)에서 추세가 휘면, 단순히 위의 모델에서 예측변수를 포함하는 \(x=t\)와 \(c=\tau\)로 바꾸어 적을 수 있습니다. \[\begin{align*} x_{1,t} & = t \\ x_{2,t} = (t-\tau)_+ &= \left\{ \begin{array}{ll} 0 & t < \tau\\ (t-\tau) & t \ge \tau \end{array}\right. \end{align*}\] \(x_{1,t}\)과 \(x_{2,t}\)의 관련 계수가 \(\beta_1\)와 \(\beta_2\)이면, \(\beta_1\)는 \(\tau\) 시점 전의 추세의 기울기, \(\tau\) 시점 이후의 선의 기울기는 \(\beta_1+\beta_2\)입니다. 형태 \((t-\tau)_+\)에 추가적인 변수를 넣어 관계식에 구부러지는 부분을 더 넣을 수 있습니다. 여기에서 \(\tau\)는 “매듭”이나 선이 휘는 시점입니다.

예제: 보스턴 마라톤 우승 기록

그림 5.20의 위쪽 패널은 1897년 시작한 이래 보스턴 마라톤 우승 기록(분 단위)을 나타냅니다. 연도에 따라 우승 기록이 나아지고 있기 때문에 시계열은 일반적인 하향 추세를 나타냅니다. 아래쪽 패널은 선형 추세를 데이터에 맞춘 것으로부터 얻은 잔차를 나타냅니다. 그래프에서 선형 추세가 잡아내지 못한 분명한 비선형 패턴이 보입니다. 어느 정도의 시간에 따라 감소하는 이분산성(heteroscedasticity)도 있습니다.

Figure 5.20: 선형 추세로 보스턴 마라톤 우승 기록을 맞추는 것은 적절하지 않습니다.

\(y\) 변수를 변환하여, 맞춘 모델이 \[ \log y_t=\beta_0+\beta_1 t +\varepsilon_t \] 가 되도록 지수 추세를 (로그-선형 회귀와 같은) 데이터에 맞출 수 있습니다. 이것 역시 이분산성(heteroscedasticity)을 설명합니다. 맞춘 지수 추세와 예상값을 그림 5.21에 나타냈습니다. 지수 추세가 선형 추세보다 데이터에 훨씬 더 잘 맞는 것 같지는 않지만, 우승 기록이 고정된 선형적 비율보다는 어떤 감소율을 가지고 감소할 것이라는 것을 나타내는 결과를 좀 더 알아볼 수 있는 형태로 냅니다.

기록 그래프에서 다른 시기 3개가 보입니다. 처음부터 약 1940년까지는 우승 기록의 변동성이 큽니다. 전반적으로 감소합니다만 1920년대에 유의미하게 증가합니다. 1940년 이후에는 시간에 따라 거의 선형적으로 감소하다가 1980년대 이후 평평해지고 그리고 표본의 끝까지 거의 상승하는 모습으로 나타납니다. 여기에서 조심해야 할 부분을 하나 강조하겠습니다. 데이터의 적당한 매듭(knot)을 주관적으로 찾는 것은 모델의 예측 성능에 해를 끼칠 수 있는 과도한 맞춤으로 이어질 수 있으니 주의 깊게 수행해야 합니다.

h <- 10

fit.lin <- tslm(marathon ~ trend)
fcasts.lin <- forecast(fit.lin, h = h)

fit.exp <- tslm(marathon ~ trend, lambda = 0)
fcasts.exp <- forecast(fit.exp, h = h)

t <- time(marathon)
t.break1 <- 1940
t.break2 <- 1980

tb1 <- ts(pmax(0, t - t.break1), start = 1897)
tb2 <- ts(pmax(0, t - t.break2), start = 1897)

fit.pw <- tslm(marathon ~ t + tb1 + tb2)

t.new   <- t[length(t)] + seq(h)
tb1.new <- tb1[length(tb1)] + seq(h)
tb2.new <- tb2[length(tb2)] + seq(h)

newdata <- cbind(t=t.new, tb1=tb1.new, tb2=tb2.new) %>%
  as.data.frame()
fcasts.pw <- forecast(fit.pw,newdata = newdata)

fit.spline <- tslm(marathon ~ t + I(t^2) + I(t^3) +
  I(tb1^3) + I(tb2^3))
fcasts.spl <- forecast(fit.spline,newdata = newdata)

autoplot(marathon) +
  autolayer(fitted(fit.lin), series = "선형") +
  autolayer(fitted(fit.exp), series = "지수") +
  autolayer(fitted(fit.pw), series = "조각별") +
  autolayer(fitted(fit.spline), series = "3차 스플라인") +
  autolayer(fcasts.pw, series="조각별") +
  autolayer(fcasts.lin, series="선형", PI=FALSE) +
  autolayer(fcasts.exp, series="지수", PI=FALSE) +
  autolayer(fcasts.spl, series="3차 스플라인", PI=FALSE) +
  xlab("연도") +  ylab("우승 기록 (분)") +
  ggtitle("보스턴 마라톤") +
  guides(colour=guide_legend(title = " "))

Figure 5.21: 보스턴 마라톤 우승 기록에 대한 선형, 지수적, 조각별 선형, 3차 스플라인 추세로부터 얻은 예측값을 추정하기.

위의 그림 5.21 은 적합선과 선형, 지수, 조각별 선형, 3차 스플라인 추세로 얻은 예측값을 나타냅니다. 조각별 선형으로 얻은 것이 가장 좋은 것 같습니다. 3차 스플라인이 과거 데이터에 가장 잘 맞지만 예측값은 좋지 않습니다.

3차 스플라인이 끝 부분에서 선형이 되도록 하는 몇 가지 제한조건이 더 들어가는 대안이 있는데 “자연적 3차 평활 스플라인(natural cubic smoothing splines)”이라고 합니다. 보통 이렇게 변형한 것이 과거 데이터를 맞추는 정도의 변화 없이 종종 더 나은 예측값을 냅니다. 그림 5.22 에서 3차 스플라인 예측값을 얻기 위해 splinef() 함수를 사용했습니다. 그림 5.21 에서 사용한 것보다 더 많은 매듭(knot)을 사용합니다만, 과대적합(over-fitting)을 방지하기 위해 계수에 제한조건이 걸리고 양 끝에서 곡선이 선형이 됩니다. 여기에는 매듭 선택 방식이 주관적이지 않게 된다는 장점도 있습니다. 이분산성(heteroscedasticity)을 다루기 위해 로그 변환도 사용했습니다.

marathon %>%
  splinef(lambda=0) %>%
  autoplot() +
    xlab("연도")

Figure 5.22: 자연적 3차 평활 스플라인을 마라톤 데이터에 적용한 것. 관측된 데이터의 끝에서 추세의 선형 추정값이 예측값입니다.

그림 5.23에 나타낸 잔차에 이분산성(heteroscedasticity)이 조금 남아있지만 모델이 추세를 잘 잡아낸 것을 나타냅니다. 예측값과 관련된 넓은 예측 구간은 과거 우승 기록에서 나타난 변동성을 반영하는 조각별 선형(piecewise linear) 추세를 이룹니다.

marathon %>%
  splinef(lambda=0) %>%
  checkresiduals()

Figure 5.23: 3차 스플라인 추세로부터 얻은 잔차.