Chapter 2. OLS (with R)

 

Chapter 2 Linear Regression by OLS and MLE | Abracadabra Stats - Companion Book (bookdown.org) 이 사이트의 코드를 참고하여 글을 작성 하였습니다. 이 사이트와는 다르게 R에있는 32개의 관측치와 11개의 numeric variables 로 이루어진 mtcars 데이터를 사용했습니다.

---

이번 Chapter에서는 Chapter 1. 에서 언급한 오차제곱합을, 함수를 이용한 오차제곱합, Nelder-Mead algorithm을 이용하여 구한 오차제곱합을 비교해보도록 한다.

1. 패키지 불러오기

#install.packages("ggplot")
#install.packages("dplyr")
#install.packages("plotly")
#install.packages("optimization")

library(dplyr)
library(ggplot2)
library(plotly)
library(optimization)

 

2. 임의의 회귀직선 정하기

 먼저 설명변수 1개만 쓸 것이기 때문에 mtcars에 있는 종속변수(mpg)와 가장 상관성이 있는 변수만 이용한다.

data("mtcars")

head(mtcars)

select_Var=NA
for(i in 1:length(mtcars)-1) {
  select_Var[i] =cor(mtcars[,1],mtcars[,i+1])
}

names(mtcars[which.max(abs(select_Var))+1])
print(abs(select_Var)[5])

#> names(mtcars[which.max(abs(select_Var))+1])
#[1] "wt"

#> print(select_Var[5])
#[1] -0.867659

설명변수 자동차 중량(wt)와 종속변수 자동차 연비(mpg)의 상관성이 -0.867659로 가장높은 상관성을 가지고 있다.

먼저 이 두변수를 이용하여 산점도를 그려본다.

ggplot(mtcars, aes(x=wt,y=mpg))+geom_point()+theme_light()

 

연비(mpg)와 중량(wt) scatterplot

이 그래프를 보면 두 변수간에 음의상관성이 있는것을 알 수 있다. 또, 어림잡어 생각해 봤을 때, 기울기가 -8이고 절편이 45정도 되는 회귀직선을 만족할거 같이 보인다. 회귀직선을 \(beta_0, \beta_1\)을 최적화 하기전에 임의로 정한 회귀 직선을 그려보았다.

 

ggplot(mtcars, aes(x=wt,y=mpg)) + geom_point()+xlab("wt") +
  geom_abline(slope =-8, intercept = 45, colour = "blue") +
  geom_abline(slope = 0, intercept = mean(mtcars$mpg),colour="red")+
  theme_light()

임의로 추정한 회귀직선

회귀직선을 임의로 정해서 그런지 위 그림 오른쪽 3개의 data point는 잘 적합해 보이지 않는다. 따라서 이를 직접 확인하기 위해서 data point 들과 회귀직선간의 선을 그려본다.

 

attach(mtcars)

y1 = (wt*-8)+45
y1_resid = y1-mpg
y2 = (wt*0) + mean(mpg)
y2_resid = y2-mpg

data = cbind.data.frame(wt,mpg,y1,y1_resid,y2,y2_resid)
data

ggplot(data, aes(wt,mpg))+ geom_point() + xlab("wt") +
  geom_abline(slope = -8, intercept =45, colour = "blue") +
  geom_segment(xend = wt, yend = y1, linetype = "dotted") +
  theme_light()

그림을 보면 기울기가 더 작아져야 한다는 것을 유추해 볼 수있다.

 

3. 오차제곱합 비교

$$  \sum_i^n = (y_i- \beta_0 - \beta_1 x_1)^2 $$

 

내가 임의로만든 회귀직선과 기울기가 0일 때의 회귀직선의 오차제곱합을 비교해보자.

 

sum((y1_resid)^2)

#[1] 509.5833

sum((y2_resid)^2)

#[1] 1126.047

 내가 만든 회귀직선이 더 작은것을 확인 할 수 있다. 하지만 나는 데이터를 보고 직접 예측하여 만들었기 때문에 당연히 더 좋을 수 밖에 없다. 이제 기울기가 0인 회귀직선을 가지고 Nelder-Mead algorithm 을 이용하여 오차제곱합을 최적화 시켜보도록 하자.

 

4. Nelder-Mead algorithm

 

---

  Nelder and Mead method 설명 (tistory.com) Nelder Mead method는 이 사이트에 잘 설명이 되어있으니 참고해 보도록한다.

---

Optimization은 내가 최적화시키려는 함수를 먼저 정의하고(Objective function, Loss function, cost function) 이에 해당하는 제약을 주고싶을때는 Constraint Functions을 이용하여 제약을하게 되는데, 이 내용을 일반화하면 다음과 같다.

$$ min_{x\in D} \qquad f(x) $$

$$ subject\ to \ g_i(x) \le 0, i = 1,...m $$

$$ \qquad  h_j(x) = 0 ,j =1,...r $$

 

여기서 \(f(x)\)는 Objective function 이고 subject to 뒤에있는 함수들이 Constraint Functions 이다.

 

이 글의 목적은 오차제곱합 \(\sum_i^n (y_i - \beta_0 -\beta_1 x_1)^2 \)을 최소로하는 \(\beta_0, \beta_1\)을 찾는 것이다. 즉, 기울기\(\beta_1\), 절편\(\beta_0\)를 찾는 것이다.

 

최적화를 하는과정에서 여러 알고리즘을 이용하는데 여기서는 Nelder-Mead algorithm을 이용해보도록 하겠다.

 

leastsq <- function(x){ 
  sum((y - x[1] - x[2]*wt)^2) #sum(y - B_0 - B_1*x_i)^2
}

먼저 최적화하려는 Objective function을 정의해준다.

 

Output<-optim_nm(fun = leastsq, k = 2,trace=TRUE)
plot(Output)

위 사진은 Nelder-Mead algorithm을 수행하는 과정이다 여기서 x1이 절편, x2가 기울기에 해당된다. 이 과정을 더욱더 잘 확인하기 위해서 3d plot으로 그려봤다.

 

a<-Output$trace
a<-as.data.frame(a)

plot_ly(a, x = a$x_1, y = a$x_2, z = a$function_value)

 

Nelder-Mead Algorithm 과정

 

Output$par

# 37.285429 -5.344383 

오차제곱합은 \(\beta_0\)가 37.285429, \(\beta_1\)이 -5.344383 일때 최소가 되는것을 알 수 있다. 추정된 모수들의 값으로 그래프를 그려 확인해보자.

 

ggplot(mtcars, aes(x=wt,y=mpg)) + geom_point()+xlab("wt") +
  geom_abline(slope =Output$par[2], intercept = Output$par[1], colour = "blue") +
  theme_light()

그래프상으로 임의로 적합한 회귀직선보다 잘 적합된것을 확인할 수 있다.

 

y3 = (wt*Output$par[2])+Output$par[1]
y3_resid = y3-mpg

sum((y3_resid)^2)

# [1] 278.3219

수치상으로도 오차제곱합이 매우많이 줄어 들었다.

 

5. lm함수를 이용한 오차제곱합

R에서는 회귀분석을 할때 알아서 model parmeter을 찾아주고 적합해주는 lm함수가 있다.

이 함수를 이용하고 Nelder-Mead algorithm의 오차제곱합 값과 비교해 보자.

m1 = lm(mpg ~ wt, data = mtcars)
m1$coefficients

#(Intercept)          wt 
#  37.285126   -5.344472 

Chapter 1. 단순회귀분석 (python) (tistory.com)에서 언급했던 방법으로 직접 계산을 해보면

 

X = model.matrix(mpg ~ wt)
solve(t(X)%*% X)%*% t(X)%*%mpg

#                 [,1]
#(Intercept) 37.285126
#wt          -5.344472

두 값이 같다는 것을 확인 할 수있다. lm를 이용하여 적합시켰을때와 Nelder-Mead algorithm을 이용하였을때의 \(\beta_0\), \(\beta_1\)의 값이 차이가 많이나지 않는 것을 알 수 있다.

 

lm함수를 이용하여 적합시켰을 때의 그래프를 확인해보자.

ggplot(mtcars, aes(x=wt,y=mpg)) + geom_point()+xlab("wt") +
  geom_abline(slope =m1$coefficients[2], intercept = m1$coefficients[1], colour = "blue") +
  theme_light()

 

그래프에서도 알고리즘과 차이가 없다는 것을 확인 할 수 있다.

y4 =(wt*m1$coefficients[2])+m1$coefficients[1]
y4_resid = y4-mpg

rbind(sum((y3_resid)^2),sum((y4_resid)^2))
#y3_resid : nelder-mead method
#y4_resid : OLS

#          [,1]
#[1,] 278.3219
#[2,] 278.3219

위 코드의 결과를 보면 오차제곱합의 차이가 없다는 것을 알 수 있다.

 

6. 결론

이 포스트에서는 오차제곱합을 Nelder-Mead 방법과, 최소제곱합(lm 함수 이용)을 이용하여 각각 결과를 확인하고 비교를 해봤다. 생각한 것과 같이 최적화를 할때 큰제약이 없는 Nelder-Mead method를 이용하면 최소제곱합방법과 비교했을 때 차이점이 없는 것을 확인할 수 있었다. 만약 다른 최적화 알고리즘(BFGS, Conjugate Gradient 등)을 이용하면 많은 제약과 초기값 설정에 있어서 어려움이 있기 때문에, 최소제곱합방법이 더 정확할 것으로 보인다.