다양한 R 스타일 가이드


좋은 코딩 스타일은 정확한 구두점을 찍는 것과 비슷합니다. 구두점이 없더라도 의미를 알 수는 있지만 구두점은 의미 전달을 더욱 쉽게 하고 글을 읽기 쉽게 만듭니다. R 에는 아래와 같이 유명한 스타일이 존재합니다.


Bioconductor’s coding standards (https://bioconductor.org/developers/how-to/coding-style/)

- Hadley Wickham’s style guide (http://stat405.had.co.nz/r-style.html )

- Google’s R style guide (http://google-styleguide.googlecode.com/ svn/trunk/google-r-style.html )

- Colin Gillespie’s R style guide (http://csgillespie.wordpress.com/2010/ 11/23/r-style-guide/)

The State of Naming Conventions in R (by Rasmus Bååth)


일관적으로 코딩하는 것은 1) 사람들과 소통하기 쉽도록하고, 2) 스스로도 읽기 쉬운 코드를 작성할 수 있습니다. 코드를 작성한지 오랜 시간이 지나고, 다시 그 코드를 보았을 때, 의미를 알기 힘든 경우가 많습니다. 하지만 최대한 비슷한 스타일로 코딩을 한다면, 코드 이해에 드는 부담을 줄일 수 있습니다. 또한 많은 사람들이 비슷한 코딩 규약을 사용한다면 코딩을 통한 소통이 편리할 것입니다.


해들리 위컴의 R 코딩 스타일 가이드


본 포스팅에서는 tidyverse (tidyr, dplyr, ggplot2 등) 패키지로 유명한 해들리 위컴의 R 스타일 가이드를 정리해보도록 하겠습니다. 해들리 위컴의 R 스타일은 Google 의 R 스타일 가이드 를 약간 변형한 스타일입니다.  


이름짓기

- 파일 이름

1. 파일 이름은 의미가 있어야합니다.

# Good
fit-models.R
utility-functions.R

# Bad
foo.r
stuff.r

2. 만약 파일이 순서대로 실행되어야 한다면, 이와 같이 숫자로 prefix 를 붙이는 것도 좋습니다. 

0-download.R
1-parse.R
2-explore.R

- 객체 이름 

“There are only two hard things in Computer Science: cache invalidation and naming things.”
Phil Karlton

변수명과 함수이름은 소문자를 사용하며, 단어를 분리할 때는 underscore(_) 를 사용합니다. 일반적으로 변수는 명사, 함수는 동사를 사용하는 것이 좋습니다. 이름은 최대한 간결하며 의미를 잘 전달할 수 있습니다. (이것은 때론 쉽지 않습니다.)

# Good
day_one
day_1

# Bad
first_day_of_the_month
DayOne
dayone
djm1

c, mean 과 같은 이미 사용되고 있는 객체에는 새로운 변수나 함수를 할당하는 것이 좋지 않습니다. 

# Bad
T <- FALSE
c <- 10
mean <- function(x) sum(x)

문법

Spacing

infix operators (=, +, -, <-, etc.) 에는 space 를 양옆으로 위치시킵니다. = 에서도 마찬가지입니다. comma 뒤는 한 칸 띄우며, comma 전에는 띄우지 않습니다. (이는 영어, 한글에서도 마찬가지입니다.)

# Good
average <- mean(feet / 12 + inches, na.rm = TRUE)

# Bad
average<-mean(feet/12+inches,na.rm=TRUE)

하지만 : 나 ::는 spacing을 하지 않습니다. 

# Good
x <- 1:10
base::get

# Bad
x <- 1 : 10
base :: get

if, for 문 등에서 ( 전에 한 칸을 띄웁니다. 하지만 함수를 호출할 때는 띄우지 않습니다. 

# Good
if (debug) do(x)
plot(x, y)

# Bad
if(debug)do(x)
plot (x, y)

아래처럼 정렬을 위해 space 를 두 번 이상하는것을 허용합니다. 

list(
  total   = a + b + c, 
  mean  = (a + b + c) / n
)

중괄호나 대괄호에는 comma 가 있는 경우가 아니라면 space 를 넣지 않습니다. 

# Good
if (debug) do(x)
diamonds[5, ]

# Bad
if ( debug ) do(x)  # No spaces around debug
x[1,]   # Needs a space after the comma
x[1 ,]  # Space goes after comma not before

중괄호

중괄호의 시작 ( { ) 이 있는 경우, 항상 개행을 해야합니다. 중괄호의 끝( } ) 이 있는 경우, if 문의 else 가 뒤따라오지 않는다면 그 줄에는 } 만 오도록 합니다. 물론 중괄호 안에서의 indentation 은 기본입니다. 

# Good

if (y < 0 && debug) {
  message("Y is negative")
}

if (y == 0) {
  log(x)
} else {
  y ^ x
}

# Bad

if (y < 0 && debug)
message("Y is negative")

if (y == 0) {
  log(x)
} 
else {
  y ^ x
}

매우 짧은 if 문의 경우 한 줄에 작성하는 것을 허용합니다. 

if (y < 0 && debug) message("Y is negative")

한 줄의 길이

한 줄의 길이를 80 자가 넘지 않도록 하는 것이 좋습니다. 이는 일반적인 크기의 글꼴로 프린트를 했을 때 용지에 잘 맞도록 하기 위함입니다. 

Indentation

Indentation 을 할 때, space 만 사용하는 것이 좋습니다. 다만 function 을 정의할 때, parameter 의 수가 많아 줄이 길어진다면 함수의 정의가 시작되는 부분에 맞게 indentation 을 합니다. 

long_function_name <- function(a = "a long argument", 
                               b = "another argument",
                               c = "another long argument") {
  # As usual code is indented by two spaces.
}

할당

변수 할당에는 = 를 사용하지 않고 <- 를 사용합니다. 

# Good
x <- 5
# Bad
x = 5

주석 달기

코드에 주석을 다는 것이 좋습니다. 주석은 무엇을 (what) 이 아닌 왜(why) 에 대한 내용을 설명해야 합니다. what 에 대한 내용은 코드 자체로 충분히 이해할 수 있도록 하는것이 좋습니다. 아래처럼 주석을 입력한 줄에 - 또는 = 을 사용해 읽기 쉽도록 분리하는 것도 좋습니다. 

# Load data ---------------------------

# Plot data ---------------------------


출처 - http://adv-r.had.co.nz/Style.html




PEP 은 무엇인가?


PEP 8 -- Style Guide for Python Code

PEP:8
Title:Style Guide for Python Code
Author:Guido van Rossum <guido at python.org>, Barry Warsaw <barry at python.org>, Nick Coghlan <ncoghlan at gmail.com>
Status:Active
Type:Process
Created:05-Jul-2001
Post-History:05-Jul-2001, 01-Aug-2013

PEP8은  Python Enhancement Proposals (PEP) 의 8 번째 내용으로, 더 나은 파이썬 코드를 작성하기 위한 하나의 코딩 규약 (Style Guide for Python Code)입니다. (https://www.python.org/dev/peps/pep-0008/)


Python PEP8 code convention cheat sheet


아래 파이썬 코드는 PEP8 style 을 활용한 코딩 예입니다. 일종의 cheat sheet로 이 코드만 알면 PEP8 에서 지향하는 파이썬 코딩 스타일의 대부분을 알 수 있습니다!


내용을 간단하게 요약하면 아래와 같습니다.


1. 모듈과 패키지 이름은 짧고 lower_case_with_underscore 이다. 

2. 다른 모듈에서 import 할 때 보일 필요가 없는 함수, 변수는 변수명 앞에 _를 추가한다. 

3. 상수는 A_CONSTANT 

4. 함수는 naming_convention

5. 클래스는 NamingConvention

6. 한 라인의 길이가 79 길이가 넘지 않도록 한다. (이것은 대략적인 A4 용지 사이즈이다.)


#! /usr/bin/env python # -*- coding: utf-8 -*- """This module's docstring summary line. This is a multi-line docstring. Paragraphs are separated with blank lines. Lines conform to 79-column limit. Module and packages names should be short, lower_case_with_underscores. Notice that this in not PEP8-cheatsheet.py Seriously, use flake8. Atom.io with https://atom.io/packages/linter-flake8 is awesome! See http://www.python.org/dev/peps/pep-0008/ for more PEP-8 details """ import os # STD lib imports first import sys # alphabetical import some_third_party_lib # 3rd party stuff next import some_third_party_other_lib # alphabetical import local_stuff # local stuff last import more_local_stuff import dont_import_two, modules_in_one_line # IMPORTANT! from pyflakes_cannot_handle import * # and there are other reasons it should be avoided # noqa # Using # noqa in the line above avoids flake8 warnings about line length! _a_global_var = 2 # so it won't get imported by 'from foo import *' _b_global_var = 3 A_CONSTANT = 'ugh.' # 2 empty lines between top-level funcs + classes def naming_convention(): """Write docstrings for ALL public classes, funcs and methods.     Functions use snake_case.     """ if x == 4: # x is blue <== USEFUL 1-liner comment (2 spaces before #) x, y = y, x # inverse x and y <== USELESS COMMENT (1 space after #) c = (a + b) * (a - b) # operator spacing should improve readability. dict['key'] = dict[0] = {'x': 2, 'cat': 'not a dog'} class NamingConvention(object): """First line of a docstring is short and next to the quotes.     Class and exception names are CapWords.     Closing quotes are on their own line     """ a = 2 b = 4 _internal_variable = 3 class_ = 'foo' # trailing underscore to avoid conflict with builtin # this will trigger name mangling to further discourage use from outside # this is also very useful if you intend your class to be subclassed, and # the children might also use the same var name for something else; e.g. # for simple variables like 'a' above. Name mangling will ensure that # *your* a and the children's a will not collide. __internal_var = 4 # NEVER use double leading and trailing underscores for your own names __nooooooodontdoit__ = 0 # don't call anything (because some fonts are hard to distiguish): l = 1 O = 2 I = 3 # some examples of how to wrap code to conform to 79-columns limit: def __init__(self, width, height, color='black', emphasis=None, highlight=0): if width == 0 and height == 0 and \ color == 'red' and emphasis == 'strong' or \ highlight > 100: raise ValueError('sorry, you lose') if width == 0 and height == 0 and (color == 'red' or emphasis is None): raise ValueError("I don't think so -- values are %s, %s" % (width, height)) Blob.__init__(self, width, height, color, emphasis, highlight) # empty lines within method to enhance readability; no set rule short_foo_dict = {'loooooooooooooooooooong_element_name': 'cat', 'other_element': 'dog'} long_foo_dict_with_many_elements = { 'foo': 'cat', 'bar': 'dog' } # 1 empty line between in-class def'ns def foo_method(self, x, y=None): """Method and function names are lower_case_with_underscores.         Always use self as first arg.         """ pass @classmethod def bar(cls): """Use cls!""" pass # a 79-char ruler: # 34567891123456789212345678931234567894123456789512345678961234567897123456789 """ Common naming convention names: snake_case MACRO_CASE camelCase CapWords """ # Newline at end of file


출처 - https://gist.github.com/RichardBronosky/454964087739a449da04



R - dplyr 을 통한 데이터 변형과 장점


본 포스팅에서는 tidyverse의 핵심 구성원이라고 할 수 있는 dplyr 을 통해 데이터를 변형하는 기본적인 방법을 basic R과 비교하여 설명하겠습니다


dplyr 의 5가지 핵심 함수

  • -값을 기준으로 선택하라 (filter)
  • -행을 재정렬하라 (arrange)
  • -이름으로 변수를 선택하라 (select)
  • -기존 변수들의 함수로 새로운 변수를 생성하라 (mutate)
  • -많은 값을 하나의 요약값으로 합쳐라 (summarize)

사용할 데이터셋은 뉴욕시에서 2013년에 출발한 336,776 개의 모든 항공편이 포함된 데이터 (nycflights13 패키지의 flights 데이터셋)

library(nycflights13)
library(tidyverse)
nrow(flights)
## [1] 336776

Filter

  • -1월 1일만 고르기
jan <- flights[flights$month == 1 & flights$day == 1, ]
nrow(jan)
## [1] 842
jan <- filter(flights, month == 1, day == 1)
nrow(jan)
## [1] 842
  • -출발 지연 시간이 120 미만인 항공편 고르기
  • -여기서 basic R과 dplyr 의 차이가 드러나는데, basic R 의 경우 arr_delay 변수가 NA 인 경우, 모든 column이 NA인 row를 포함한 dataframe을 반환한다.
  • -basic R의 경우 filter 를 구현할 때, NA 인 경우도 고려해야하므로, dplyr이 더 간단하게 원하는 목적을 달성할 수 있다고 볼 수 있다. 
'under_two <- flights[(flights$arr_delay < 120), ]
nrow(under_two) # arr_delay가 NA인 row가 모든 column이 NA인 row로 추가된다. .
## [1] 326576
sum(is.na(under_two$arr_delay))
## [1] 9430
sum(is.na(under_two$dep_delay))
## [1] 9430
under_two <- filter(flights, arr_delay < 120) 
nrow(under_two) 
## [1] 317146
sum(is.na(under_two$arr_delay))
## [1] 0
sum(is.na(under_two$dep_delay)) 
## [1] 0
# R 기본문법으로 dplyr의 filter를 구현하려면 이렇게 함. 
under_two <- flights[(flights$arr_delay < 120) & !is.na(flights$arr_delay), ]
nrow(under_two) 
## [1] 317146
sum(is.na(under_two$arr_delay))
## [1] 0
sum(is.na(under_two$dep_delay))
## [1] 0
  • -출발 혹은 도착에서 2시간 이상 지연되지 않은 항공편을 모두 찾기
  • -이 때, 두 조건절이 모두 NA 인 경우, basic R 에서는 모든 컬럼이 NA 인 행을 dataframe에 추가하여 반환한다. 
under_two <- flights[!((flights$arr_delay > 120) | (flights$dep_delay > 120)), ]
nrow(under_two) # 둘다 NA인 행 9430개의 행이, 포함되어있음
## [1] 325354
sum(is.na(under_two$arr_delay))
## [1] 9304
sum(is.na(under_two$dep_delay))
## [1] 9304
under_two <- filter(flights, !(arr_delay > 120 | dep_delay > 120)) 
nrow(under_two) # 둘다 NA인 행 제외. 이게 정상적인 결과.
## [1] 316050
sum(is.na(under_two$arr_delay))
## [1] 0
sum(is.na(under_two$dep_delay)) 
## [1] 0
  • -출발 시간, 도착 시간 둘 중 하나가 2시간 이상 지연된 항공편을 모두 찾기
over_two <- flights[((flights$arr_delay > 120) | (flights$dep_delay > 120)), ]
nrow(over_two) # 둘중 하나가 NA 이면서, delay > 120 인 행 포함, 근데 둘 다 NA인 경우 모든 컬럼이 NA인 행을 포함함 (9403 개)
## [1] 20726
sum(is.na(over_two$arr_delay)) # 9430 -> 둘중 하나가 NA인 수
## [1] 9430
sum(is.na(over_two$dep_delay)) #9304 -> 둘다 NA인 수 
## [1] 9304
over_two <- filter(flights, ((arr_delay > 120) | (dep_delay > 120))) # 둘중 하나가 NA 이면서, delay > 120인 행 포함
nrow(over_two) # 정상적인 결과 
## [1] 11422
sum(is.na(over_two$arr_delay)) 
## [1] 126
sum(is.na(over_two$dep_delay)) 
## [1] 0
  • -일반 데이터 프레임의 문제점. 조건이 NA인 경우 모든 컬럼이 NA인 행이 출력된다. 
  • -dplyr 의 filter 는 이러한 문제 없이 원하는 결과를 출력해준다. 
flights[NA,]
## # A tibble: 336,776 x 19
##     year month   day dep_time sched_dep_time dep_delay arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>
##  1    NA    NA    NA       NA             NA        NA       NA
##  2    NA    NA    NA       NA             NA        NA       NA
##  3    NA    NA    NA       NA             NA        NA       NA
##  4    NA    NA    NA       NA             NA        NA       NA
##  5    NA    NA    NA       NA             NA        NA       NA
##  6    NA    NA    NA       NA             NA        NA       NA
##  7    NA    NA    NA       NA             NA        NA       NA
##  8    NA    NA    NA       NA             NA        NA       NA
##  9    NA    NA    NA       NA             NA        NA       NA
## 10    NA    NA    NA       NA             NA        NA       NA
## # ... with 336,766 more rows, and 12 more variables: sched_arr_time <int>,
## #   arr_delay <dbl>, carrier <chr>, flight <int>, tailnum <chr>,
## #   origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>, hour <dbl>,
## #   minute <dbl>, time_hour <dttm>

Arrange

  • -날짜별로 오름차순 정렬하기
  • -정렬의 경우, 마찬가지로 dplyr 이 훨씬 직관적이다.
arranged <- flights[order(flights$year, flights$month, flights$day),]
head(arranged)
## # A tibble: 6 x 19
##    year month   day dep_time sched_dep_time dep_delay arr_time
##   <int> <int> <int>    <int>          <int>     <dbl>    <int>
## 1  2013     1     1      517            515         2      830
## 2  2013     1     1      533            529         4      850
## 3  2013     1     1      542            540         2      923
## 4  2013     1     1      544            545        -1     1004
## 5  2013     1     1      554            600        -6      812
## 6  2013     1     1      554            558        -4      740
## # ... with 12 more variables: sched_arr_time <int>, arr_delay <dbl>,
## #   carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>,
## #   time_hour <dttm>
arranged <- arrange(flights, year, month, day)
head(arranged)
## # A tibble: 6 x 19
##    year month   day dep_time sched_dep_time dep_delay arr_time
##   <int> <int> <int>    <int>          <int>     <dbl>    <int>
## 1  2013     1     1      517            515         2      830
## 2  2013     1     1      533            529         4      850
## 3  2013     1     1      542            540         2      923
## 4  2013     1     1      544            545        -1     1004
## 5  2013     1     1      554            600        -6      812
## 6  2013     1     1      554            558        -4      740
## # ... with 12 more variables: sched_arr_time <int>, arr_delay <dbl>,
## #   carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>,
## #   time_hour <dttm>
  • -내림차순 정렬하기
desc <- flights[order(-flights$arr_delay),]
head(desc)
## # A tibble: 6 x 19
##    year month   day dep_time sched_dep_time dep_delay arr_time
##   <int> <int> <int>    <int>          <int>     <dbl>    <int>
## 1  2013     1     9      641            900      1301     1242
## 2  2013     6    15     1432           1935      1137     1607
## 3  2013     1    10     1121           1635      1126     1239
## 4  2013     9    20     1139           1845      1014     1457
## 5  2013     7    22      845           1600      1005     1044
## 6  2013     4    10     1100           1900       960     1342
## # ... with 12 more variables: sched_arr_time <int>, arr_delay <dbl>,
## #   carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>,
## #   time_hour <dttm>
desc <- arrange(flights, desc(arr_delay))
head(desc)
## # A tibble: 6 x 19
##    year month   day dep_time sched_dep_time dep_delay arr_time
##   <int> <int> <int>    <int>          <int>     <dbl>    <int>
## 1  2013     1     9      641            900      1301     1242
## 2  2013     6    15     1432           1935      1137     1607
## 3  2013     1    10     1121           1635      1126     1239
## 4  2013     9    20     1139           1845      1014     1457
## 5  2013     7    22      845           1600      1005     1044
## 6  2013     4    10     1100           1900       960     1342
## # ... with 12 more variables: sched_arr_time <int>, arr_delay <dbl>,
## #   carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>,
## #   time_hour <dttm>

Select

  • -select 함수는 열을 선택하는 함수이다.
selected <- flights[, c("year", "month", "day")]
head(selected)
## # A tibble: 6 x 3
##    year month   day
##   <int> <int> <int>
## 1  2013     1     1
## 2  2013     1     1
## 3  2013     1     1
## 4  2013     1     1
## 5  2013     1     1
## 6  2013     1     1
selected <- select(flights, year, month, day)
head(selected)
## # A tibble: 6 x 3
##    year month   day
##   <int> <int> <int>
## 1  2013     1     1
## 2  2013     1     1
## 3  2013     1     1
## 4  2013     1     1
## 5  2013     1     1
## 6  2013     1     1

Mutate

  • -mutate는 새로운 열을 추가하는 함수이다.
  • -dplyr은 select를 할 때 end_with function 을 통해 해당 문자로 끝나는 컬럼을 선택할 수도 있고, : 를 통해 두 문자의 사이에 있는 컬럼을 선택할 수도 있다. 
  • -따라서 열 선택에 있어 더욱 많은 기능을 짧은 코드로 구현할 수 있으며, 새로운 변수를 정의할 때도, 한 문장에 정의할 수 있어 코드의 길이가 짧아지고 가독성이 좋아진다. 
flights_sml <- select(flights, year:day, ends_with("delay"), distance, air_time)
flights_sml_new <- mutate(flights_sml, gain=arr_delay - dep_delay, speed = distance / air_time * 60)
nrow(flights_sml_new)
## [1] 336776
  • -새 변수만을 남기고 싶다면 transmute를 사용한다. 
flights_transmute <- transmute(flights, gain = arr_delay - dep_delay,
                               hours = air_time / 60, 
                               gain_per_hour = gain / hours)
nrow(flights_transmute)
## [1] 336776

Summarize

  • -데이터 프레임을 하나의 행으로 축약한다. 
# dplyr을 사용한 방법
summarize(flights, delay = mean(dep_delay, na.rm = TRUE), delay_sd = sd(dep_delay, na.rm=TRUE))
## # A tibble: 1 x 2
##   delay delay_sd
##   <dbl>    <dbl>
## 1  12.6     40.2
  • -summarize 함수는 group_by (dplyr 에서 제공하는 함수) 와 보통 함께 사용하는 경우가 많다. 이것은 dplyr 의 이전 버전이라고 할 수 있는 plyr에서도 잘 구현되어 있었는데, dplyr 이 되면서 조금 더 이해하기 쉽게 문법이 바뀌었다. 
by_day <- group_by(flights, year, month, day)
summarize(by_day, delay = mean(dep_delay, na.rm = TRUE))
## # A tibble: 365 x 4
## # Groups:   year, month [12]
##     year month   day delay
##    <int> <int> <int> <dbl>
##  1  2013     1     1 11.5 
##  2  2013     1     2 13.9 
##  3  2013     1     3 11.0 
##  4  2013     1     4  8.95
##  5  2013     1     5  5.73
##  6  2013     1     6  7.15
##  7  2013     1     7  5.42
##  8  2013     1     8  2.55
##  9  2013     1     9  2.28
## 10  2013     1    10  2.84
## # ... with 355 more rows

파이프로 여러 작업 결합하기

  • -파이프 명령어인 %>% 를 통해 dplyr 을 통한 데이터 변형 작업을 더욱 직관적이고, 가독성있게 수행할 수 있다. 
  • -각 위치에 대해 거리와 평균 지연 사이의 관계를 탐색하고 싶은 경우, 파이프를 사용하지 않는 방법은 아래와 같다. 
by_dest <- group_by(flights, dest)
delay <- summarize(by_dest, count = n(),
                   dist = mean(distance, na.rm = TRUE),
                   delay = mean(arr_delay, na.rm = TRUE))
delay <- filter(delay, count > 20, dest != "HNL")
ggplot(data = delay, mapping = aes(x = dist, y = delay)) + 
  geom_point(aes(size = count), alpha = 1/3) + 
  geom_smooth(se=FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

  • -파이프를 사용하면? 
  • -아래와 같이 임시 변수를 따로 놓지 않으면서 단 한 번의 logic 으로 똑같은 목적을 달성할 수 있다! 
  • -이러한 파이프 방법은 일련의 데이터 변형 작업을 의미 단위로 묶어 코드가 모듈화가 되기 때문에 유저가 코드를 쉽게 구현하고 이해할 수 있도록 도와준다. 
  • -다만 파이프가 너무 길어지면 유용한 임시변수를 두는 방법이 효율적일 수 있다. 
delay <- flights %>% 
  group_by(dest) %>%
  summarize(count = n(),
            dist = mean(distance, na.rm=TRUE),
            delay = mean(arr_delay, na.rm=TRUE)
            ) %>% 
  filter(count > 20, dest != "HNL")
ggplot(data = delay, mapping = aes(x = dist, y = delay)) + 
  geom_point(aes(size = count), alpha = 1/3) + 
  geom_smooth(se=FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'


윈도우에서 R 업데이트 하기 


아래 방법으로 최신 버전의 R을 설치할 수 있습니다. 


install.packages("installr")

library(installr)

updateR()


이전 버전의 패키지 가져오기


1. 이전 R 버전 폴더의 모든 라이브러리들을 새로운 R 버전 폴더로 copy and paste 합니다.

저는 R 3.4.* 버전에서 R 3.5.3. 버전으로 업데이트하였는데, 경로는 아래와 같습니다. 


C:\Users\Username\Documents\R\win-library\3.4

C:\Users\Username\Documents\R\win-library\3.5


3.4 폴더에 있는 모든 라이브러리 폴더들을 3.5 폴더로 복사 붙여넣기하면 됩니다. 하지만 이렇게만 하면 패키지 버전이 호환되지 않기 때문에 제대로 로드가 안될 수도 있습니다. 따라서 3.5 폴더에 있는 라이브러리들을 업데이트 시켜줘야합니다.


2. 패키지 업데이트


update.packages(checkBuilt=TRUE)


를 입력하면 팝업 박스가 나오면서 yes no 를 선택할 수 있는데 yes 를 계속해서 클릭해줍니다. 


패키지 ‘ucminf’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘triebeard’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘ordinal’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘magic’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘lpSolve’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘RcppProgress’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘desc’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘xopen’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다


...... 


그러면 이런식으로 패키지를 최신 버전으로 업데이트를 알아서 해줍니다. 


3. 패키지가 모두 제대로 설치가 되었는지 확인합니다. 


version

packageStatus()


이제 최신 버전의 R을 이용할 수 있습니다. 

R - 다중회귀분석 (Multiple Linear Regression)


R 을 통해 다중회귀분석을 수행하는 기본적인 절차를 포스팅합니다. 


라이브러리 임포트


datarium 의 marketing 이라고 하는 데이터셋을 이용하여 다중회귀분석을 수행한다.

library(devtools)
# install_github("kassambara/datarium")

library(tidyverse)
library(datarium)

data("marketing", package = "datarium")
head(marketing, 4)



탐색적 자료분석 - 데이터 분포 살펴보기


가장 먼저 할 일은 데이터의 분포를 살펴보는 것이다. 

ggplot(data = marketing, aes(x = newspaper, y = sales)) +
  geom_point(size=5)




위와 같이 하나하나 산점도를 볼 수 있지만, 아래처럼 모든 변수들의 분포 (Variation)와 두 변수의 관계 (Covariation) 한 번에 볼 수 있도록 하는 라이브러리들이 많이 존재한다. 예를 들어, PerformanceAnalytics 라이브러리를 이용하면 아래와같이 한 번에 분포와 상관계수까지 구할 수 있다. 


library("PerformanceAnalytics")
chart.Correlation(marketing, histogram=TRUE, pch=19)



모델 만들기


sales = b0 + b1youtube + b2facebook + b3*newspaper 


분포를 봤을 때, youtube, facebook, newspaper 가 sales와 선형관계를 갖고 있다고 가정하고 위와 같이 다중 회귀 모델을 만들 수 있으며, 아래 R 코드를 통해 모델을 적합시킬 수 있다. 

model <- lm(sales ~ youtube + facebook + newspaper, data = marketing) summary(model)

Call:
lm(formula = sales ~ youtube + facebook + newspaper, data = marketing)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.5932  -1.0690   0.2902   1.4272   3.3951 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.526667   0.374290   9.422   <2e-16 ***
youtube      0.045765   0.001395  32.809   <2e-16 ***
facebook     0.188530   0.008611  21.893   <2e-16 ***
newspaper   -0.001037   0.005871  -0.177     0.86    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.023 on 196 degrees of freedom
Multiple R-squared:  0.8972,	Adjusted R-squared:  0.8956 
F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16


결과 해석


Regression 결과를 해석하는데, 가장 먼저할 것은 F-statistics 를 확인하는 것이다. 가장 밑에 F-statistics 를 보면, p-value가 매우 작아, 최소 한 변수가 유의하다는 것을 알 수 있다. 어떤 변수가 유의한지를 확인하기 위해, 아래 coefficient table 을 보면 된다. 이 때, t-value와 p-value 는 귀무가설 평균=0 로 놓고 구한 값으로, 귀무가설로부터 얼마나 벗어난지에 관한 지표이다. 


X의 estimate 값은 다른 변수가 고정되어있을 때, X가 한 단위 변화했을 때, Y의 평균적인 변화량이다. 예를 들어, youtube의 estimate 는 0.04 인데, 이는 youtube 가 1 변화했을 때, sales의 평균적인 변화량이다. newpaper 의 경우, p-value > 0.05 이므로 이러한 estimate 의 값이 유의하지 않다는 것을 알 수 있다. 


t value 와 p-value 는 무엇인가?


$$ H_0 : \beta = 0 $$

$$ t = \frac{\hat{\beta}}{se(\hat{\beta})} $$ 


이 때 위의 t 값은 자유도가 n-2 인 t 분포를 따른다. (이렇게 되는 이유)

p-value는 이렇게 계수 추정치가 t 분포를 따른다는 사실을 이용해서 구한 값이다. 


결정계수와 F-stat 


다음과 같이 정의해보자.


$$ SSE = \sum^{n}_{i=1}{(y_i - \hat{y})^2} $$

$$ SSR = \sum^{n}_{i=1}{(\bar{y} - \hat{y_i})^2} $$

$$ SST = \sum^{n}_{i=1}{(y_i - \bar{y})^2} $$


이 때, 각각의 통계량에 대한 자유도는 아래와 같다. 


df(SSE) : n-p

df(SSR) : p-1

df(SST) : n-1


이 때, p는 intercept를 포함한 모델의 총 파라미터의 갯수이며, n은 샘플 수이다. 


$$ MSE = \sum^{n}_{i=1}{(y_i - \hat{y})^2} / (n-p) $$

$$ MSR = \sum^{n}_{i=1}{(\bar{y} - \hat{y_i})^2} / (p-1) $$

$$ MST = \sum^{n}_{i=1}{(y_i - \bar{y})^2} / (n-1) $$


귀무가설을 아래와 같이 설정하자


$$ H_0 : \beta_1 = \beta_2 = \beta_3 = ... = \beta_{p-1} $$ 


이 때, 


$$ \frac{MSE}{MSR} \sim F(n-p, p-1) $$


이를 통해 F-stat 과 p-value를 구할 수 있다. 이 값은 (explained variance) / (unexplained variance) 가 F-stat 을 따른다는 사실을 이용한 검정 방법이라고 할 수 있다. 왜 이렇게 되는지에 대해서는 이 링크를 참고.



결정계수 (Coefficient of determination, Rsquare)


결정계수는 아래와 같이 정의되며, 전체 분산 중 모델에의해 설명되는 분산의 양을  의미한다. 


$$ R^2 = \frac{SSR}{SST} $$


adjusted 결정계수를 사용하기도 하는데, 이는 아래와 같이 정의된다. 


$$ \bar{R^2} = 1-(1-R^2) \frac {n-1}{n-p} $$ 


적합된 모델의 계수

summary(model)$coefficient


새로운 모델을 아래와 같이 구할 수 있다. sales = 3.5 + 0.045youtube + 0.187facebook.


새로운 모델

model  <- lm(sales ~ youtube + facebook, data = marketing)
summary(model)

Call:
lm(formula = sales ~ youtube + facebook, data = marketing)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.5572  -1.0502   0.2906   1.4049   3.3994 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.50532    0.35339   9.919   <2e-16 ***
youtube      0.04575    0.00139  32.909   <2e-16 ***
facebook     0.18799    0.00804  23.382   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.018 on 197 degrees of freedom
Multiple R-squared:  0.8972,	Adjusted R-squared:  0.8962 
F-statistic: 859.6 on 2 and 197 DF,  p-value: < 2.2e-16


신뢰구간

confint(model)

2.5 %97.5 %
(Intercept)2.808411594.20222820
youtube0.043012920.04849671
facebook0.172138770.20384969


위와 같이 계수의 신뢰구간을 출력할 수도 있다. 이는 t 분포를 이용한다. 


3일차


일정 : 하카타역 seattle's best coffee (아침식사) - 8:50분 하카타역 쿠루쿠루 버스 - 유후인 (긴린코 호수, 유노츠보 거리) - 유후인 길거 음식 (점심식사) - 쿠로가와 - 온천욕 - 우오가시 스시 (이른 저녁식사) - 캐널시티 - 이치란 라멘 (늦은 저녁 식사)


후쿠오카 여행 3일차에는 "쿠루쿠루 버스" 를 이용하여 유후인과 쿠로가와를 다녀왔습니다. 쿠루쿠루 버스는 주요 관광지를 순회하는 버스인데, 한 명에 6만원 정도에 이용했던 것 같습니다. 쿠루쿠루 버스는 꽤 괜찮았습니다. 산큐패스를 사더라도, 미리 버스 시간 예약하고 발권 받고, 위치 확인하고 등등 귀찮은 일이 많은데, 쿠루쿠루 버스는 정해준 시간/장소에 집결지에 가기만 하면 알아서 관광지를 데려다 주니까 여행 일정 짜는데 시간이 부족하신 분들은 정말 강추합니다! 다만 한 여행지를 느긋하게 보기는 힘들다는 단점은 있었습니다. 


유후인 긴린코 호수 & 유노츠보 거리

 



제가 여행했던 날 (2019년 3월 1일) 은 정말 날씨가 좋았습니다. 덕분에 긴린코 호수에서 호수 배경의 산들이 잘 나오게 사진이 찍혀서 만족스러웠습니다. 긴린코 호수는 상상했던것 보다 작았지만 그래도 한적하게 산책하기 정말 좋았습니다. 좋은 풍경에 공기도 맑아, 걷기만해도 좋은 기분이었습니다.  









길거리 음식으로 당고를 사먹었습니다. 오른쪽에 있는 것이 기본 당고이고, 왼쪽은 김이랑 마요네즈가 토핑된 것인데 둘다 떡이 부드럽고 맛있었어요.



유후인에 가면 꼭 먹는다는 금상고로케인데, 정말 명성대로 맛있었습니다.







어떤 도시를 여행을 하고 나서 얼마 후, 다시 또 그곳에 가고 싶은 기분이 든다면, 자신과 잘 맞는 좋은 여행지라 생각합니다. 제게 유후인은 그런 여행지 중 하나였던 것 같아요. 점심을 굶고, 길거리 음식만 먹으면서 둘러보았음에도 버스 승차 시간 때문에 느긋하게 보지 못한 게 너무 아쉽더라구요. 유후인은 다음에는 꼭 느긋하게 보고 싶고, 꼭 1박을 해보고 싶은 그런 여행지였습니다. 화려한 볼 거리는 없지만 왠지 모르게 정감이 가는 그런 마을이었습니다. 



쿠로가와 온천마을


구마모토에 위치한 쿠로가와는 유후인에서 버스로 50분 정도 걸려 도착했습니다. 쿠로카와에서는 2시간 정도의 관광 시간이 주어졌는데, 이 시간 동안 온천을 먼저 하고 남은 시간에는 마을 구경을 하고 버스에 승차하기로 했습니다. 









쿠로가와 유모토소 료칸 온천


쿠로가와에 있는 여러 온천 중 어떤 온천에 들어갈지 미리 정해오지 않아서 즉흥적으로 유모토소 온천에 들어갔습니다. 대중탕을 할지 가족탕을 할지 고민하다가 사적으로 이용할 수 있는 가족탕을 이용하였습니다. 가족탕 중 대나무 형태로 된 욕조를 선택했고, 50분 이용에 가격은 2만원 정도 됐던 것 같습니다.













나중에 안 사실이지만, 유모토소 온천은 쿠로카와 내에서도 오래된 것으로 유명한 온천이라고 하더라구요. 그래서 그런지 가족탕 내부 시설이 다소 낡아있었기 때문에 깔끔한 느낌은 아니었습니다. 전날 갔던 후쿠오카 근교의 나카가와 세이류 온천이 상대적으로 훨씬 깔끔했던것 같네요. 다만 이 온천의 장점은 대나무 숲을 바라보며 온천욕을 즐길 수 있다는 것이었습니다. 근데 문제는, 왼쪽에 사람이 지나다니는 다리가 위치해있어 문을 활짝 열고 온천을 즐기지는 못했습니다. ㅎㅎ 유모토소 온천을 즐기고 나오는데 몸에서 쇠 냄새 같은 것이 나길래, 욕조가 녹슬었나 했는데 알고보니 유모토소 온천은 온천수에 철이 많이 함유된으로 유명한 온천이라고 하더라구요. 이런 온천을 함철천이라고 부르고, 드물면서도 몸에도 좋은 온천이라고 하던데 일반적인 온천이 아니라 색다른 온천을 경험을 한 것 같아 좋았습니다. 





유모토소 온천 근처 슈크림 빵집에서 슈크림 빵을 사먹었습니다. 사람들이 되게 많이 먹고 있어 궁금해서 먹어봤는데 맛있긴 맛있었습니다. 온천하고 나와서 먹으면 기분 좋을듯한 그런 달콤한 맛이었습니다.






쿠로가와는 도심에서 떨어진 시골인데다가 산 속에 위치해 있어 도시에만 살았던 저로서는 신선한 느낌을 받았습니다. 또 관광객이 그렇게 많지 않아 일본스러운 분위기를 많이 느낀 곳이기도 합니다. 쿠로카와에 사는 주민 분들은 무슨 생각과 감정을 갖고 사는지가 정말 궁금해졌습니다. 평소 살아가는 세계어서 벗어나 새로운, 좋은 세계를 경험한 것 같아 정말 기분이 좋았습니다. 


돌아오는 버스에서는 나른해져 계속 잠만 잤습니다. 잠이 깨려고 할 때쯤 하카타 역에 도착했습니다. 


캐널시티 인형뽑기


하카타역 지하에 위치한 하카타 1번가에서 우오가시 스시를 먹고, 무엇을 하지 고민하다가 캐널시티 게임센터에 갔습니다. 평소에도 인형뽑기를 즐겨하는데, 일본 인형뽑기방은 어떨지 궁금했습니다. 가보고 딱 든 생각은 인형들이 각맞춰서 줄 서 있는 느낌이라고 할까요. ㅎㅎ 그리고 인형이 뽑히면 바로 그자리를 직원이 와서 채워놓아서 기계 안의 인형들의 형태를 항상 같게 유지시키는 것 같습니다. 





저는 이 곳에서 블루투스 이어폰을 득템했습니다. (한 판에 100엔이었고 약 2천 500엔 정도 썼습니다..) 뽑고나서 정말 기뻤는데, 문제는 이어폰이 한쪽만 들어있는 이어폰이었습니다. 이왕 만들거면 양쪽 귀에 꽂을 수 있도록 하지 왜 한쪽 이어폰만 상품으로 만들었는지 참.. 그래도 기념으로 뽑은걸로 생각하고 있습니다.



후쿠오카 이치란 라멘 본점


지난번 후쿠오카 여행 때 먹어보지 못해 이번에는 꼭 먹어보려고 계획했던 이치란 라멘입니다. 밤 9시쯤 간거같은데, 줄이 꽤 있었습니다. 자판기에 돈을 넣은 후, 원하는 음식의 티켓을 뽑고 들어가 앞에다 내놓으면 음식을 주는 시스템이었습니다. 티켓을 뽑고 마음대로 안에 들어가면 안되고 앞에서 기다리면 직원이 안내해줍니다. 





이치란 라멘은 정말 정말 맛있었습니다. 지금까지 알던 라멘맛이랑 다르게 국물 맛이 깊으면서도 느끼하지 않은 맛이었습니다. 차슈가 굉장히 맛있으니 꼭 차슈 추가를 하시고, 약간 짠 편이라 공기밥과 같이 드시는걸 추천합니다! 지난번 후쿠오카 여행 때 아무 음식점이나 들어가서 라멘을 먹었다가 실패한 경험이 있는데, 왜 이치란 라멘을 한국 사람들이 많이 가는지 알 것 같았습니다. 하카타역에서 이른 저녁을 먹고 왔음에도 불구하고, 남김 없이 깨끗하게 다 먹었습니다. 


이렇게 후쿠오카 여행 3일차 포스팅을 마치겠습니다 :)


2019 후쿠오카 여행기


2019-02-27 ~ 03-02 까지 일본 후쿠오카 여행을 다녀왔습니다! 항공권을 늦게 끊어 조금 비싸게 구입 (왕복 27만원..) 해서 전체 예산은 항공권 포함 80만원 정도로 잡았습니다. 후쿠오카는 온천과 맛있는 음식들이 유명하기 때문에 그 위주로 여행을 하였습니다. 


1일차


일정 : 8시경 일본 도착 - 오오야마 모츠나베 - 편의점 - 숙소




오오야마 모츠나베


오오야마 모츠나베는 맛을 3가지 중에 하나를 선택할 수 있었는데, 된장맛을 선택했는데 맛있었습니다. 약간 짠 편이라 흰 쌀밥이 꼭 필요한 맛입니다. 밥 시켜 드시길 추천 드립니다. 다만 이 집은 한국인 입맛에 잘 맞는 모츠나베 집으로 유명한 편이라 한국 사람들이 엄청 많았습니다. 화장실 가려고 가게 내부를 둘러보니 거의 70% 정도가 한국인이었습니다. 후쿠오카의 대부분의 음식점들이 그렇지만 한국어 메뉴판이 구비되어있고, 종업원들도 서비스가 좋았습니다.



저녁 비행기였고, 밤에 비까지와서 첫 날 일정은 모츠나베로 마무리했습니다. 호텔에 돌아오면서 편의점에서 파르페와 푸딩을 사와서 먹었습니다. (호텔은 하카타 엑셀 도큐 호텔로 선택했는데, 나카스, 캐널시티, 텐진이 모두 가까워서 정말 위치적으로 좋았습니다!)


2일차


일정 : 텐진 역 안 도토루 카페 (아침식사) - 오호리 공원 - JR 하카타 시티 우마야 (점심 식사) - 나카가와 세이류 - 텐진 우시부치 (저녁식사) - 야쿠인 - 돈키호테 - 나카스강 - 숙소



텐진역 도토루


2일차에는 텐진역사 내부에 있는 도토루 카페에서 아침식사를 하였습니다. 가격은 4~5천원 정도 했던것 같습니다. 일본 여행하면서 아침 식사로 뭐먹을지 고민하시는 분들께, 도토루 카페에서 아침 식사 해보는 것을 추천드립니다. 가격대비 너무 훌륭한 아침식사였습니다. 도토루는 일본의 흔한 카페인데도 불구하고 빵의 퀄리티가 매우 좋고 드립 커피도 깔끔하게 맛있었습니다. 





다음으로 오호리 공원을 갔는데, 아쉽게도 비가 조금 와서 많이 둘러보지는 못하고 한 30분 둘러보다 나왔습니다.. 다음에 후쿠오카에 가서 날씨 좋으면 꼭 다시 들러보고 싶네요. 






아뮤 프라자에 있는 포켓몬 스토어에 들렀습니다.





JR 하카타시티 우마야


점심으로 무엇을 먹을지 하카타여 주변을 배회하다가 JR 하카타 시티에 있는 우마야에서 점심식사를 하였습니다. 닭요리를 전문으로 하는 음식점인것 같은데, 이 곳을 선택한 이유는 토리히츠마부시를 먹어보기 위해서였습니다. 토리히츠마부시는 나고야식 장어덮밥인 히츠마부시의 닭 버전으로, 훈제 닭고기 덮밥 정도로 보시면 될 것 같습니다. 다만, 차를 부어 말아먹는 것이 특이한 점입니다. 덮밥에 차를 말아먹는게 생소하긴 했지만, 맛은 정말 맛있었습니다. 







나카가와 세이류 온천


후쿠오카 근교 나카가와 세이류 온천에 다녀왔습니다. 나카가와 세이류 온천은 왕복 2시간 정도로 온천을 즐길 수 있다는 점이 좋았고, 시설이 깔끔하다는 평이 많아서 선택하였습니다. 가족탕으로 히노키탕을 선택해서 이용했는데, 물도 맑고 깔끔하고, 간단하게 온천을 즐기고 나오기 좋았습니다. 


제 경우에는 동생이 일본에 있어 전화로 예약하여 이용하였습니다. 이용하고 싶은 시간에 웹 페이지에서는 예약이 완료되었다고 표시되었는데도 불구하고, 전화로 하니까 원하는 시간에 예약이 되었습니다. 전화로도 예약을 받으니 가능하다면 전화로 예약하는 것이 확실하고 좋은 방법일듯 합니다. 원하는 탕과 시간, 이름, 전화번호 정도만 말하면 되니 일본를 하지 못하더라도 어렵지 않을 것 같습니다. 예약을 하였으면 예약 시간에 가서 예약자명을 대면 이용 방법 안내를 해줍니다. 물론 대중탕 이용하실분들은 예약이 필요 없습니다. 


또한 나카가와 세이류는 무료 셔틀 버스를 운행하는데 (https://m.blog.naver.com/PostView.nhn?blogId=obscura1&logNo=221248318706&proxyReferer=https%3A%2F%2Fwww.google.com%2F) 이 포스트를 참고하면 하카타에서 무료 셔틀버스를 이용하는 방법을 알 수 있습니다! 셔틀버스 시간표는 자주 바뀌는듯하니 최신 시간표를 참고하세요.


나카가와 세이류 인터넷 예약 : https://nakagawaseiryu.resv.jp/reserve/calendar.php?x=1552105410

무료 셔틀버스 시간표 : http://www.nakagawaseiryu.jp/pdf/bus_time.pdf








텐진 우시부치


일본 야끼니쿠를 꼭 먹어보고 싶었는데, 이번에 그 소망을 이루었습니다. 우시부치는 텐진에 위치한 야끼니쿠 전문점인데, 분위기가 굉장히 좋고, 세련된 음식점이었습니다. 먼저 모듬세트 2인분 (약 2만 6천원)을 시켜 먹었는데 신기하게도, 원래 알고 있던 소고기 맛이랑 달랐습니다. 고기를 얇게 썰어서 그런지 몰라도 같은 소고기인데 훨씬 부드러웠습니다. 각 부위별로 이름표가 나와서 무엇인지 확인하면서 먹는 재미도 있었습니다. 모듬을 다 먹고, 양이 조금 부족해서 맛있게 먹은 부챗살 (1만 5천원)을 한 번 더 시켜서 먹었습니다. 고기 5점에 1만 5천원이었는데, 부챗살이 비싼 부위는 아닌것을 감안하면 싼 건 아닌 것 같습니다. 가성비를 중요시 여기신다면 모듬으로 시켜먹는것을 추천드립니다. 이 집의 야끼니쿠가 이번 여행에서 가장 맛있었던 것 같습니다.







텐진에서 숙소로 걸어 돌아오면서 나카스 야타이를 잠깐 들렀습니다.




 오늘 하루 이렇게 먹고도 돌아오면서 닭껍질 꼬치와 맥주를 사와서 먹었습니다. 


Logistic regression의 maximum likelihood estimation 직관적 이해



Logistic regression의 Random component 



$$ Y \sim B(n, \pi) $$ 


Logistic regression은 반응변수 y가 위와 같은 이항 분포를 따른다고 가정한다.  



Model 


$$ \pi = \frac{1}{1+e^{-{(X\beta+\alpha)}}} $$ 


이 때, X를 통해 모수 pi 를 위와 같이 추정하는 것이다. 또는 일반적으로 아래와같이 표현되기도 한다. 


$$ logit(\pi) = X\beta + \alpha $$ 


Logistic regression에서 하고자하는건 결국, Random component에 주어진 pi를 추정하는 것인데, 이것이 0-1 사이의 값을 갖기 때문에 logit 변환을 수행하는 것이다. logit link function은 0~1사이의 값을 음의 무한대에서 양의 무한대로 바꾸어준다. 


Random component가 이항 분포이므로 likelihood는 아래와 같이 계산된다.



이곳에 Model을 집어넣음으로써, beta와 alpha 를 maximum likelihood 추정법으로 계산할 수 있게 된다. 근데 여기서 n과 y는 n번 시도했을 때 성공횟수 y를 나타낸다. 근데, 실제 각각의 데이터를 보면 시행횟수 n=1이며, y=0 or 1이다. 따라서 위 likelihood 는 베르누이 likelihood 와 같게된다. 



결국, 위의 likelihood를 최대화하면 계수 추정값을 얻을 수 있는데, log likelihood를 최대화하는 방식으로 쉽게 계수를 추정할 수 있다. 


$$ \sum_{i=1}^{N} y_i log(\pi_i) + (1-y_i) log (1-\pi_i) $$


Logistic regression의 log likelihood이다. 이는 negative binary cross entropy 이다. 결국, log likelihood를 최대화 하는 것은 negative binary cross entropy를 최대화 하는 것과 같고 이는 binary cross entropy를 최소화화는 것과 같다. 

고급통계 - Net reclassification improvements


Net reclassification index란 예측 모형에서 새로운 모델이 과거 모델보다 얼마나 얼마나 더 나아졌는지를 측정하는 지표이다. AUC 값을 기준으로 이를 판단하는게 과거에는 일반적이었으나, 그 한계점이 지적되외서 NRI가 등장하게 되었다. Framingham Heart Study에서 HDL 관련 연구를 수행할 때, AUC의 한계를 볼 수 있었는데, HDL을 예측 모형에 추가했을 때, AUC값이 크게 증가하지 않는 것이었다. 하지만 AUC 측면에서가 아니라, 모형의 계수를 통해 HDL의 effect를 파악할 때, HDL은 심장질환에 중요한 예측 인자였다. 


AUC의 한계는 다음과 같이 제시되었다. 1) 임상적인 의미와 연결시키기 힘들다. 2) 작은 변화는 잡아내기가 힘들다.


하지만 NRI는 이전 모델과 비교해 새로운 모델에서 얼마나 많은 사람들이 제대로 재분류되었는가를 파악함으로써 모델의 개선을 측정한다. 

NRI example table
EventTest 1Total, splitTotal
Non-eventAbnormalNormal
Test 2Abnormal1842228
246
Normal26872
85664
Total, split201030
106070
Total3070100

출처 - https://en.wikipedia.org/wiki/Net_reclassification_improvement


만약 환자 30명, 비환자 70명에 대해서 환자-비환자 여부를 판단하는 모델 Test1, Test2를 만들었다고 하자. Test 2가 Test1에 비해 얼마나 우월한지를 판단하고자 한다. 위 테이블에서 여기서 검은색 글씨는 Test 1, Test 2 둘 다 맞춘거고, 하얀색 글씨는 Test1, Test 2 둘 다 틀린 경우, 초록색 글씨는 Test 2에서만 맞추고, Test 1에서는 틀린 것, 빨간색 글씨는 Test 1에서 맞추고, Test 2 에서 틑린 것이다. 


이 때, NRI의 계산은 매우 간단한데, Test 2를 통해서 얻은 이득이 얼마나 더 많은지를 계산해주면 된다. 이것은 Test 2를 도입했을 때, 추가적으로 맞춘것-오히려 틀린것 을 해줌으로써 계산되며, Event, Non-event 두 그룹의 대상자에 대해서 계산한 다음에 더함으로써 최종 NRI를 계산하게 된다. 


즉, Event가 발생한 사람들에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_e = (4-2)/30 = 0.067 $$ 

다음으로 Non-event에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_{ne} = (8-4)/70 = 0.057 $$ 

따라서 NRI는 아래와 같이 계산된다.

$$ NRI = NRI_e + NRI_{ne}= 0.124 $$ 


Gradient Boosting Algorithm의 직관적인 이해


실패를 통해 성공을 발전시켜라. 낙담과 실패는 성공으로 가는 가장 확실한 두 개의 디딤돌이다. -데일 카네기


Gradient Boosting Algorithm (GBM)은 회귀분석 또는 분류 분석을 수행할 수 있는 예측모형이며 예측모형의 앙상블 방법론 중 부스팅 계열에 속하는 알고리즘입니다. Gradient Boosting Algorithm은 Tabular format 데이터 (엑셀형태와 같이 X-Y Grid로 되어있는 데이터)에 대한 예측에서 엄청난 성능을 보여주고, 머신러닝 알고리즘 중에서도 가장 예측 성능이 높다고 알려진 알고리즘입니다. 그렇기 때문에 Gradient Boosting Algorithm을 구현한 패키지들이 많습니다. LightGBM, CatBoost, XGBoost 같은 파이썬 패키지들이 모두 Gradient Boosting Algorithm을 구현한 패키지들입니다. GBM은 계산량이 상당히 많이 필요한 알고리즘이기 때문에, 이를 하드웨어 효율적으로 구현하는 것이 필요한데, 위 패키지들은 모두 GBM을 효율적으로 구현하려고한 패키지들이라고 볼 수 있습니다.


Boosting 이란?


주요 앙상블 알고리즘은 bagging과 boosting으로 나눌 수 있고, Gradient boosting은 boosting 계열의 앙상블 알고리즘입니다. Gradient Boosting 알고리즘은 Gradient 를 이용하여 Boosting하는 알고리즘입니다. Gradient boosting 외의 boosting 계열의 알고리즘의 예로는 Adaptive boosting을 들 수 있습니다. 


먼저 Boosting의 개념부터 살펴봅시다. Boosting이란 약한 분류기를 결합하여 강한 분류기를 만드는 과정입니다. 분류기 A, B, C 가 있고, 각각의 0.3 정도의 accuracy를 보여준다고 합시다. A, B, C를 결합하여 더 높은 정확도, 예를 들어 0.7 정도의 accuracy를 얻는 게 앙상블 알고리즘의 기본 원리입니다. Boosting은 이 과정을 순차적으로 실행합니다. A 분류기를 만든 후, 그 정보를 바탕으로 B 분류기를 만들고, 다시 그 정보를 바탕으로 C 분류기를 만듭니다. 그리고 최종적으로 만들어진 분류기들을 모두 결합하여 최종 모델을 만드는 것이 Boosting의 원리입니다. 



GBM의 직관적인 이해


GBM을 이해하는 가장 쉬운 방법은 Residual fitting으로 이해하는 것입니다. 아주 간단한 모델 A를 통해 y를 예측하고 남은 잔차 (residual)을 다시 B라는 모델을 통해 예측하고 A+B 모델을 통해 y를 예측한다면 A보다 나은 B 모델을 만들 수 있게 되죠. 이러한 방법을 계속하면 잔차는 계속해서 줄어들게되고, training set을 잘 설명하는 예측 모형을 만들 수 있게 됩니다. 하지만 이러한 방식은 bias는 상당히 줄일 수 있어도, 과적합이 일어날 수도 있다는 단점이 있습니다. 따라서 실제로 GBM을 사용할 때는 sampling, penalizing 등의 regularization 테크닉을 이용하여 더 advanced 된 모델을 이용하는 것이 보편적입니다. 하지만 본 포스팅의 목적은 GBM 에 대하여 직관적인 이해를 해보는 것이기 때문에 이 부분은 다음에 기회가 있을 때 정리해보도록 하겠습니다. 

위 그림을 보시면 tree 1을 통해 예측하고 남은 잔차를 tree2를 통해 예측하고, 이를 반복함으로서 점점 잔차를 줄여나가는 것을 볼 수 있습니다. 이 때, 각각의 모델 tree1,2,3 을약한 분류기 (weak learner), 이를 결합한 분류기를 강한 분류기 (strong learner)라고도 합니다. 보통 약한 분류기로는 간단한 의사결정나무 (decision tree)를 많이 사용합니다. 이를 Gradient boosting tree라고도 하는데, 구현한 대표적인 라이브러리로 XGboost를 들 수 있습니다. XGBoost는 python, R로 이용할 수 있습니다. 



Gradient란 무엇인가?


그렇다면 residual fitting과 gradient가 무슨 관계인지 궁금하실텐데요. residual은 loss function을 squared error로 설정하였을 때, negative gradient 입니다. 따라서 residual에 fitting해서 다음 모델을 순차적으로 만들어 나가는 것 negative gradient를 이용해 다음 모델을 순차적으로 만들어 나가는 것으로 볼 수 있습니다. 우선 아래 식을 통해 이를 확인해봅시다. 


loss function을 아래와 같이 정의하면,


$$ j(y_i, f(x_i)) = \frac{1}{2}(y_i - f(x_i))^2 $$ 


이 때의 negative gradient는 residual이 됩니다.


residual이 $$ y_i - f(x_i) $$ 이기 때문에 negative gradient가 residual이 된다는 것을 알 수 있습니다. 


따라서 gradient boosting은 다음 모델을 만들 때, negative gradient를 이용해 만들기 때문에 gradient boosting 이라고 할 수 있습니다. 또한 residual fitting model은 gradient boosting 모델의 한 가지 종류입니다. 왜냐하면 다른 loss function을 이용해서 gradient boosting을 할 수 있기 때문입니다. 예를 들어 회귀 문제가 아닌 분류 문제라면 loss function으로 squared error 를 이용할 수 없기 때문에 새로운 loss function을 만들어야하는데, 이 경우에도 negative gradient를 이용해서 새로운 model 을 fitting하고 이를 합산해 나가시는 식으로 최종 모델을 만들게 됩니다. 


직관적으로 이것은 무엇을 뜻할까요? 왜 negative gradient를 이용해서 새로운 모델을 만드는 걸까요? negative gradient는 pseudo-residual이라고도 불리며, 이것은 어떤 데이터 포인트에서 loss function이 줄어들기 위해 f(x)가 가려고하는 방향입니다. 이 방향에 새로운 모델을 fitting해서 이것을 이전 모델과 결합하면, f(x) 는 loss function이 줄어드는 방향으로 업데이트가 되겠죠. 이것이 gradient boosting의 아이디어입니다. Gradient boosting을 gradient descent + boosting 이라고 하기도 합니다. 어쨌든, loss function을 줄이는 방향의 negative gradient를 얻고, 이를 활용해 boosting을 하는 것이기 때문에 gradient descent와 boosting이 결합된 방법이다. 라고 이해하셔도 괜찮습니다. 



위 그림은 위키피디아에서 gradient boosting에 대하여 설명한 알고리즘도입니다. 이제 위 내용을 이해할 수 있습니다. gradient boosting에서는 learning rate를 통하여 pseudo-residual에 fitting된 모델을 어느정도로 update할지에 대해서 정하게 되는데, 이 값은 위 알고리즘도의 3번처럼, loss를 최소화하는 값을 optimization 과정을 통해 얻을 수도 있겠지만, 일반적으로 0.001~0.01 정도로 매우 작은 값으로 설정하는 것이 보편적입니다. 왜냐하면 작은 값으로 설정하여야 굉장히 세밀한 분류기를 얻을 수 있습니다. learning rate가 높을 수록 빠르게 모델의 bias를 줄여나가지만, learning rate가 적으면 디테일한 부분을 놓칠 수 있다고 이해하시면 됩니다. 이와 관련해서는 gradient boosting 과정을 simulation 해볼 수 있는 사이트를 소개합니다. 이 곳에서 직접 gradient boosting 방법이 어떻게 함수를 fitting 해나가는지를 시각적으로 확인할 수 있습니다. 


참고

https://machinelearningmastery.com/gentle-introduction-gradient-boosting-algorithm-machine-learning/