[Data Mining with R] R / Simple Linear Regression - Part 1.

by Geol Choi | November 23, 2014

이번 글에서는 Linear Regression에 대한 기초 통계 이론에 대한 소개와 이에 대한 R 프로그래밍에 대해 알아보기로 한다.

Linear Regression은 간단하게 말해, 관찰된 데이터들의 변수들 간 관계를 1차원적인 Graph로 표현(이를 fitting이라고 함)하는 것이다. Linear Regression은 통계학의 역사관점에서 볼 때, 특정 변수가 다른 변수와 어떤 상관관계인지를 알아보기 위한 수단으로 발전해 왔다. 데이터를 관찰하여 이에 대한 모델을 세우고 이 모델을 통해 데이터에 대한 예측을 하고자 하는 것이 목표이며, 더 나아가 이에 대한 신뢰도를 어떻게 평가할 수 있는가가 이 이론에 대한 거의 전부라고 할 수 있다.

물론 모든 데이터의 경향이 1차원적이지 않다. 아니 1차원적이 아닌 것이 훨씬 많다. 다만, 경향성을 평가하기 위한 모델을 공부하는 데 있어 Linear Regression의 이해는 필수적이며, 향후 Non-linear Regression이라든가 Multi-variable Linear Regression으로의 이해를 하기 위한 기본 단계라고 보면 된다. 다중 변수에 대한 Linear Regression과 구분짓기 위하여 단일 변수에 대해서는 Simple Linear Regression이라고 칭하도록 한다. 그러나 이 글에서는 편의상 "Simple"을 생략하기로 한다.

Theoretical Background

다음과 같이 n개의 관찰된 데이터가 있다고 하자:

\((x_1,y_1),...,(x_n,y_n)\)

여기서 \(y_i\)는 랜덤변수 \(Y_i\)로부터 관찰된 데이터로서 기울기(slope) \(\beta_1\)와 절편(Intercept) \(\beta_0\)를 갖는 Linear Regression 모델은

\(y_i=\beta_0+\beta_1x_i+\epsilon_i\)

와 같이 표현되며, Error Term \(\epsilon_{i}\)는 평균 0, 표준편차 \(\sigma^{2}\)의 정규분포를 따른다. 즉,

\(\epsilon_{i}\sim N(0,\sigma^{2})\)

이며, 실제 관측된 데이터 \(y_{i}\)와 Linear Regression 모델을 이용하여 Fitting된 \(\hat{y}_{i}\)의 차이를 의미한다. 즉,

\(\epsilon_{i}=y_{i}-\hat{y}_{i}=y_{i}-(\beta_{0}+\beta_{1}x_{i})\)

또한 랜덤변수 \(Y_{i}\)는 다음과 같은 정규분포를 따른다.

\(Y_{i}\sim N(\beta_{0}+\beta_{1}x_{i},\sigma^{2})\)

그러면 \(\beta_{0}\)와 \(\beta_{1}\)는 어떻게 구해지는지 알아보자.

Basic Idea는 매우 단순하다. 각 변수 \(X_{i}\)에 대하여 관측된 값 \(y_i\)와 Fitting된 값 \(\hat{y}_{i}\)의 차이의 제곱의 합을 최소화하는 Line을 찾아내는 것이다. 제곱을 하는 이유는, 두 값의 상대적 차이에 따라 0보다 크거나 작을 수도 있기 때문이며, 이러한 모델을 Least Squares Model이라고 한다. 이를 수학적으로 표현하면,

\(Q=\sum_{i}^{n}\epsilon_{i}^{2}=\sum_{i}^{n}[y_{i}-(\beta_{0}+\beta_{1}x_{i})]^2\)

Q를 Linear Squares Fit이라고 한다. 이를 최소화하는 \(\beta_{0}\) 및 \(\beta_{1}\)는 Q를 각각에 대해 편미분한 값이 0이 되도록 하여 구할 수 있다. (이에 대한 자세한 내용은 Optimisation 이론을 참고하기 바란다. 사실 편미분한 값이 0인 위치에서 최소값을 갖는 것에 대해서는 Convexity 조건을 만족하는지도 확인하는 것이 맞으나 일단 이 조건에 대한 설명은 생략하도록 한다.)

\(\frac{\partial Q}{\partial \beta_{0}}=-\sum_{i}^{n}2[y_{i}-(\beta_{0}-\beta_{1}x_{i})]=0\)(*)

그리고

\(\frac{\partial Q}{\partial \beta_{0}}=\sum_{i}^{n}2x_{i}[y_{i}-(\beta_{0}+\beta_{1}x_{i})]=0\) (**)

이며, (*)과 (**)식을 각각 정리하면,

\(n \beta_{0}+(\sum_{i}^{n}x_{i})\beta_{1}=\sum_{i}^{n}y_{i}\) (***)

\((\sum_{i}^{n}x_{i})\beta_{0}+(\sum_{i}^{n}x_{i}^{2}\beta_{1})=\sum_{i}^{n}x_{i}y_{i}\) (****)

(***)과 (****)을 연립하여 풀면,

\(\hat{\beta_{1}}=\frac{n\sum_{i}^{n}x_{i}y_{i}-(\sum_{i}^{n}x_{i})(\sum_{i}^{n}y_{i})}{n\sum_{i}^{n}x_{i}^{2}-(\sum_{i}^{n}x_{i})^2}\)

\(\hat{\beta_{0}}=\frac{\sum_{i}^{n}{y_{i}}}{n} - \hat{\beta_{1}}\frac{\sum_{i}^{n}{x_{i}}}{n} = \bar{y}-\hat{\beta_{1}\bar{x}}\)

가 된다. 여기서 \(\hat{\beta_{0}}\)와 \(\hat{\beta_{1}}\)는 계산된(Estimate) 값을 의미한다. 한편,

\(S_{xx}=\sum_{i}^{n}{(x_{i}-\bar{x})^2}=\sum_{i}^{n}{x_{i}^{2}-n\bar{x}^{2}}=\sum_{i}^{n}{(\frac{\sum_{i}^{n}{x_{i}}}{n})^2}\)

및

\(S_{XY}=\sum_{i}^{n}{(x_{i}-\bar{x})(y_{i}-\bar{y})}=\sum_{i}^{n}{x_{i}y_{i}-n\bar{x}\bar{y}}=\sum_{i}^{n}{\frac{(\sum_{i}^{n}{x_{i}})(\sum_{i}^{n}{y_{i}})}{n}}\)

이므로, \(\hat{\beta_{1}}\)는 다음과 같이 다시 쓸 수 있다.

\(\hat{\beta_{1}}=\frac{S_{XY}}{S_{XX}}\)

오차에 대한 제곱합 SSE는 다음과 같이 정의된다:

\(\mathbf{SSE}=\sum_{i}^{n}{\epsilon_i^2}=\sum_{i}^{n}{(y_{i}-\hat{y}_{i})^2}=\sum_{i}^{n}{(y_{i}-(\beta_{0}+\beta_{1}x_{i}))^2}\)

이를 풀어쓰면,

\(\mathbf{SSE}=\sum_{i}^{n}{y_{i}^{2}}-\hat{\beta_{0}}\sum_{i}^{n}{}-\hat{\beta_{1}}\sum_{i}^{n}{x_{i}y_{i}}\)

이며, 오차 분산 \(\hat{\sigma}_{0}^{2}\)은

\(\hat{\sigma}_{0}^{2}=\frac{\mathbf{SSE}}{n-2}=\frac{\sum_{i}^{n}{y_{i}^{2}}-\beta_{0}\sum_{i}^{n}{y_{i}}-\hat{\beta_{1}\sum_{i}^{n}{x_{i}y_{i}}}}{n-2}\)

이다. n이 아닌 (n-2)로 나누는 이유는 Unbiased Estimator를 위한 것인데, 이에 대한 자세한 내용은 여기를 참고하기 바란다. 지금까지의 파라미터 추정을 계산하는 방법을 정리하면 다음과 같다:

\(S_{XX}=\sum_{i}^{n}{x_i^2}-{\begin{pmatrix}\sum_{i}^{n}{x_i}\end{pmatrix}^2}/{n}\)
\(S_{XY}=\sum_{i}^{n}{x_{i}y_{i}}-{(\sum_{i}^{n}{x_i})(\sum_{i}^{n}{y_i})}/{n}\)
\(\hat{\beta}_{1}=S_{XY}/S_{XX}\)
\(\hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1}\hat{x}\)
\(\mathbf{SSE}=\sum_{i}^{n}{y_i^2}-\hat{\beta}_{0}\sum_{i}^{n}{y_i}-\hat{\beta}_{1}\sum_{i}^{n}{x_i y_i}\)
\(\hat{\sigma}^2=\mathbf{SSE}/(n-2)\)

R Code

이제 이론적 배경을 토대로 Linear Regression에 대한 R 코드를 작성해 보도록 하자.

linearRegression.R:

# @Title: Linear Regression
# @Author: Geol Choi (cinema4dr12@gmail.com)
# @Purpose: To calculate intercept(beta0) and slope(beta1) for linear regression model
# and also sum of square error
linearRegression = function( x, y ) {
  sum_x = sum(x);
  sum_y = sum(y);
  sum_x2 = sum(x * x);
  sum_y2 = sum(y * y);
  sum_xy = sum(x * y);
  
  # number of data points
  n = length(x);
  
  Sxx = sum_x2 - (sum_x)^2/n;
  Sxy = sum_xy - (sum_x * sum_y)/n;
  
  # slope of linear regression model
  beta1 = Sxy/Sxx;
  
  # intercept of linear regression model
  beta0 = mean(y) - beta1 * mean(x);
  
  # residual error
  err = (sum_y2 - beta0*sum_y - beta1*sum_xy) / (n-2);
  
  return(c(beta0, beta1, err));
}




예를 들어, R의 빌트-인 데이터 중 하나인 mtcars의 mpg(Miles/Gallon)와 disp(Displacement)에 대하여 Linear Regression 모델을 계산해 보자. 우선 mpg와 disp간의 관계를 한 눈에 확인할 수 있도록 그래프로 표현하면 다음과 같다:

두 변수의 관계는 대체적으로 음의 상관관계임을 알 수 있다. 이에 대한 R 코드는 다음과 같다:

myPlot.R:

myPlot = function() {
  library(ggplot2);
  
  ggplot() + 
    coord_cartesian() +
    scale_x_continuous() +
    scale_y_continuous() +
    layer(
      data = mtcars,
      mapping = aes(x = mpg, y = disp, colour = factor(disp)),
      stat = "identity",
      geom = "point",
      size = 3
    )
}

위의 코드를 수행하기 전 주의할 점이 있다. ggplot2라는 R 패키지가 설치되어야 실행된다는 것이다. CRAN을 통해 ggplot2를 설치하려면 R 콘솔에 다음과 같이 명령을 실행한다:

> install.packages("ggplot2")

이제 두 변수 간 Linear Regression 모델을 계산해 보자:

R의 콘솔 입력창에 다음과 같이 입력한다:

> lr = linearRegression(mtcars$mpg, mtcars$disp)
> lr
[1]  580.88382  -17.42912 4470.68729




결과는 lr이라는 변수에 저장하였으며, 순서대로 \(\beta_{0}\)(Intercept), \(\beta)_{1}\)(Slope), SSE이다.

예상대로 Slope가 (-), 즉 음의 상관관계임을 알 수 있다. 이에 대한 모델은 다음과 같이 표현할 수 있다:

\(y=580.88-17.43x+\epsilon\)

그러면 이렇게 계산된 모델이 데이터 포인트들과 겹쳐서 표현되면 어떻게 나오는지 확인해 보자:

Linear Regression 잘 계산되었음을 알 수 있다. 이에 대한 코드는 다음과 같다(코드가 썩 잘 정리되지는 않은 것 같다):

myPlot.R:

myPlot = function() {
  library(ggplot2);
  
  # linear regression model
  x1 = min(mtcars$mpg);
  x2 = max(mtcars$mpg);
  y1 = 580.88382 - 17.42912*x1;
  y2 = 580.88382 - 17.42912*x2;
  vec1 = c(x1,x2);
  vec2 = c(y1,y2);
  df = data.frame(x = vec1, y = vec2);
  
  ggplot() + 
    coord_cartesian() +
    scale_x_continuous() +
    scale_y_continuous() +
    layer(
      data = mtcars,
      mapping = aes(x = mpg, y = disp, colour = factor(disp)),
      stat = "identity",
      geom = "point",
      size = 3
    ) +
    layer(
      data = df,
      mapping = aes(x = x, y = y),
      stat = "identity",
      geom = "line",
      size = 1
    )
}




그런데 R의 빌트-인 명령어로 위의 복잡한 과정을 단 한 줄에 끝낼 수 있다. 약간 허무할 수 있다고 생각할 지 모르지만 우리는 이 과정을 모두 이해했다. 모르고 사용하는 것과 알고 사용하는 것은 천지 차이이니 말이다.

R의 콘솔에 다음과 같이 입력해 보자:

> x = mtcars$mpg;
> y = mtcars$disp;
> model = lm(y ~ x);
> model

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
     580.88       -17.43 




와우~ 계산이 일치한다.

지금까지 Linear Regression 함수를 직접 만들어 사용하였는데, 이제는 R의 빌트-인 함수를 사용하여 Linear Regression의 Shiny App을 만들어 보자.

Shiny App

R Studio의 "shiny" 라이브러리를 이용하여 mtcars의 Linear Regression을 계산하는 Web 기반 실시간 어플리케이션을 만들어 보자.

우선 어플리케이션의 결과를 확인해 보자.

상기 어플리케이션의 코드는 다음과 같다. (코드에 대한 설명은 코드 내 주석으로 대신한다.)

ui.R:

# @title: Linear Regression for mtcars
# @author: Geol Choi, ph.D (cinema4dr12@gmail.com)
# @date: 2014.11.30

library(ggplot2);

shinyUI(fluidPage(
  title = "Linear Regression of mtcars",
  theme = "bootstrap.css",
  
  # Application title
  headerPanel(h2("Linear Regression by CINEMA4DR12")),
  
  # Sidebar
  sidebarPanel(
    h3("mtcars"),
    selectInput("input_x",
                label = HTML("<h4>&nbsp&nbsp&nbsp X</h4>"),
                choices = list("Miles per Gallon" = 1,
                               "Number of Cylinders" = 2,
                               "Displacement" = 3,
                               "Gross Horsepower" = 4,
                               "Rear Axle Ratio" = 5,
                               "Weight" = 6,
                               "1/4 Mile Time" = 7,
                               "V/S" = 8,
                               "Transmission" = 9,
                               "Number of Forward Gears" = 10,
                               "Number of Carburetors" = 11
                              ),
                selected = 1
    ),
    
    selectInput("input_y",
                label = HTML("<h4>&nbsp&nbsp&nbsp Y</h4>"),
                choices = list("Miles per Gallon" = 1,
                               "Number of Cylinders" = 2,
                               "Displacement" = 3,
                               "Gross Horsepower" = 4,
                               "Rear Axle Ratio" = 5,
                               "Weight" = 6,
                               "1/4 Mile Time" = 7,
                               "V/S" = 8,
                               "Transmission" = 9,
                               "Number of Forward Gears" = 10,
                               "Number of Carburetors" = 11
                ),
                selected = 1
    )
    
  ),
  
  # Show the generated 3d bivar plot
  mainPanel(
    tabsetPanel(
      tabPanel("Plot", plotOutput("plot")),
      tabPanel("Summary", verbatimTextOutput("summary"))
    )
  )
))

server.R:

# @title: Linear Regression for mtcars
# @author: Geol Choi, ph.D (cinema4dr12@gmail.com)
# @date: 2014.11.30

library(shiny);
library(ggplot2);

shinyServer(function(input, output, session) {
  Data = reactive({
    # get selection x & y from select box
    x_val = input$input_x
    y_val = input$input_y
    
    switch(x_val,
           "1" = {
             vec1 = mtcars$mpg
             name1 = "mpg"
           },
           "2" = {
             vec1 = mtcars$cyl
             name1 = "cyl"
           },
           "3" = {
             vec1 = mtcars$disp
             name1 = "disp"
           },
           "4" = {
             vec1 = mtcars$hp
             name1 = "hp"
           },
           "5" = {
             vec1 = mtcars$drat
             name1 = "drat"
           },
           "6" = {
             vec1 = mtcars$wt
             name1 = "wt"
           },
           "7" = {
             vec1 = mtcars$qsec
             name1 = "qsec"
           },
           "8" = {
             vec1 = mtcars$vs
             name1 = "vs"
           },
           "9" = {
             vec1 = mtcars$am
             name1 = "am"
           },
           "10" = {
             vec1 = mtcars$gear
             name1 = "gear"
           },
           "11" = {
             vec1 = mtcars$carb
             name1 = "carb"
           }
    )
    
    switch(y_val,
           "1" = {
             vec2 = mtcars$mpg
             name2 = "mpg"
           },
           "2" = {
             vec2 = mtcars$cyl
             name2 = "cyl"
           },
           "3" = {
             vec2 = mtcars$disp
             name2 = "disp"
           },
           "4" = {
             vec2 = mtcars$hp
             name2 = "hp"
           },
           "5" = {
             vec2 = mtcars$drat
             name2 = "drat"
           },
           "6" = {
             vec2 = mtcars$wt
             name2 = "wt"
           },
           "7" = {
             vec2 = mtcars$qsec
             name2 = "qsec"
           },
           "8" = {
             vec2 = mtcars$vs
             name2 = "vs"
           },
           "9" = {
             vec2 = mtcars$am
             name2 = "am"
           },
           "10" = {
             vec2 = mtcars$gear
             name2 = "gear"
           },
           "11" = {
             vec2 = mtcars$carb
             name2 = "carb"
           }
    )
    
    df = data.frame(tmp1 = vec1, tmp2 = vec2);
    names(df) = c(name1, name2);
    return(df);
  })
  
  output$plot = renderPlot({
    myPlot(Data())
  })
  
  output$summary = renderPrint({
    mySummary(Data())
  })
  
})


################################################
# FUNCTION : myPlot
################################################
myPlot = function(myData) {
  # extract x & y from myData
  x = myData[,1];
  y = myData[,2];
  
  # calculate linear regression from R build-in function lm(...)
  myModel = lm(y ~ x);
  coeff = coefficients(myModel);
  beta0 = as.numeric(coefficients(myModel)[1]);
  beta1 = as.numeric(coefficients(myModel)[2]);
  
  # construct a line graph from the linear regression model
  x1 = min(x);
  x2 = max(x);
  y1 = beta0 + beta1*x1;
  y2 = beta0 + beta1*x2;
  
  vec1 = c(x1,x2);
  vec2 = c(y1,y2);
  df = data.frame(x = vec1, y = vec2);
  
  name1 = names(myData)[1];
  name2 = names(myData)[2];
  aesData = aes_string(x = name1, y = name2);
  
  # plotting
  ggplot() + 
    coord_cartesian() +
    scale_x_continuous() +
    scale_y_continuous() +
    layer(
      data = mtcars,
      mapping = aesData,
      stat = "identity",
      geom = "point",
      size = 3
    ) +
    layer(
      data = df,
      mapping = aes(x = x, y = y),
      colour = "red",
      stat = "identity",
      geom = "line",
      size = 1
    )
}

################################################
# FUNCTION : mySummary
################################################
mySummary = function(myData) {
  # extract x & y from myData
  x = myData[,1];
  y = myData[,2];
  
  # calculate linear regression from R build-in function lm(...)
  myModel = lm(y ~ x);
  summary(myModel);
}

이것으로 Linear Regression Part 1.을 마치도록 하겠다. Part 2.에서는 Linear Regression 모델에 대한 파라미터 추정에 대한 이론을 다루도록 하겠다.