Hard skills/R (46)

해들리 위컴은 어떻게 수많은 R 패키지를 개발할 수 있었을까?


R 에는 수많은 패키지들이 존재한다. 이들 중 어떤 패키지를 활용해서 데이터를 분석하고 결과를 도출해야 할까? 힘드게 문법을 읽힌 패키지가 아무도 쓰지 않는 패키지 였다면? 개발자가 더 이상 유지보수를 하지 않는다면? 다른 패키지의 함수들과 호환이 되지 않는다면? 우리의 시간은 제한되어 있기에 패키지 선정에 신중해야하며 위험이 존재하기도 한다. R packages for data science 라는 슬로건을 내걸은 Tidyverse 는 R 유저들에게 고민에 빠진 R 유저들에게 그 방법을 제시한다. Tidyverse 는 같은 철학은 공유하는 데이터과학을 위한 R 패키지들의 모임이며, 함께 사용할 수 있도록 디자인되었다. 


Hadleyverse 라고도 불리는 Tidyverse 패키지들


놀라운 것은 tidyverse 에 속한 패키지들의 많은 부분이 한 명에 의해 개발되었다는 것이다. 해들리위컴은 R 패키지 그룹인 Tidyverse 에 속한 dplyr, readr, tidyr, ggplot 와 같은 패키지들의 개발자로 유명하며, 그가 R 에 기여한 공헌도는 어마어마하다. 특히 tidy data 의 개념과  Leland Wilkinson 의 The Grammar of Graphics 을 구현한 층층이 쌓아 그래프를 그리는 ggplot2 는 가장 유명하다. 


2014년 Quora 에 어떻게 해들리 위컴이 이렇게 많이 R 에 기여할 수 있었는지에 대해 질문이 올라왔는데, 해들리 위컴이 직접 답변을 달았다. 그 내용은 아래와 같다. 


1. 글쓰기 (Writing)


그는 매일 아침 60~90분 동안 글쓰기를 한다. 아침에 일어나서 가장 먼저 하는 것이 글을 쓰는 것이다. 그가 글 쓰기를 중요시하는 첫 번재 이유는 직접 쓴 글은 스스로 참고하는데 유용하기 때문이다. 예를 들어, 그는 C++ 을 항상 코딩하는 것이 아니기 때문에, 참조가 필요한데 자신이 쓴 Rcpp 에 대한 글을 참조한다고 한다. 두 번째로, 글쓰기는 지식과 Tool 의 갭을 매울 수 있고, 이는 새로운 문제를 해결하는 방법을 고안하는데 도움이 된다고 한다. 


2. 읽기 (Reading)


그는 많이 읽는다. 약 300개의 blog 를 follow 하고 있고, Twitter 나 Stack Overflow 의 R tag 가 달린 글을 자주 본다고 한다. 물론 대부분의 글은 자세히 보지는 않지만, 이러한 수많은 글에 스스로를 노출시킴으로써 기술적인 변화를 따라가는데 도움되며, 새로운 프로그래밍 언어와 다른 사람들이 R 을 통해 데이터를 다루는 방법을 꾸준히 보고 어떠한 문제가 생겼을 때, 그 '이름' 을 빠르게 인식할 수 있다. 따라서 googling 등을 통해 가능한 해결법을 찾는 것을 더욱 빠르게 한다. 


3. 청킹 (Chunking)


문맥의 전환 (Coxtext switching)  은 매우 비싼 것이다. 많은 패키지들을 한 번에 작업하면, 이도저도 아니게 된다. 따라서 그는 대부분의 시간에서는 자신의 패키지들에 관한 어떠한 이슈나 아이디어를 점점 쌓아가고, 그렇게 쌓인 것 많아졌을 때 쯤, 그 패키지를 개선하는 작업에 며칠 간의 시간을 보낸다. 


마지막으로 그는 오래전부터 R 패키지 개발에 시간을 보냈기 때문에 그러한 경험이 축적되어 새로운 패키지를 순식간에 만들어낼 수 있는 지식과 경험을 보유하고 있다고 할 수 있다. 그의 학위 논문도 데이터를 exploration 하고 modeling 하는 R 패키지에 관한 것이었다. (http://had.co.nz/thesis/) 결국 그가 지금까지 R 커뮤니티에 공헌한 것들은 이러한 수년 간의 경험과 지식의 정리, 그리고 기술적 진보를 이해하는데 들인 노력의 산물이라고 할 수 있을 것이다.  


참고

https://www.quora.com/How-is-Hadley-Wickham-able-to-contribute-so-much-to-R-particularly-in-the-form-of-packages


Apply 계열 함수 정리


 Overview

FunctionDescription
applyApply functions over array margins
byApply a function to a data frame split by factors
eapplyApply a function over values in an environment
lapplyApply a function over a list or vector
mapplyApply a function to multiple list or vector arguments
rapplyRecursively apply a function to a list
tapplyApply a function over a ragged array

출처 - https://csgillespie.github.io/efficientR/programming.html



  • 데이터를 다룰 때, 원자별, 그룹별로 함수를 적용할 경우가 많다. 
  • Apply 계열의 함수는 데이터 구조를 갖는 R object 를 input 으로 받아 원소 별 혹은 그룹별로 함수를 적용시키는 것
  • input 과 output 의 형태에 따라 여러 종류로 나뉜다.
    • apply (input : array, output : array)
    • lapply (input : list or vector, output : list)
    • sapply (input : list or vector, output : vector or array)
    • vapply (input : list or vector, output : vector or array)
    • tapply (input : list or vector and factor, output : vector or array)
    • mapply (input : list or vector, output : vector or array)

apply

x <- matrix(1:9, c(3,3))  
x 
#>      [,1] [,2] [,3]
#> [1,]    1    4    7
#> [2,]    2    5    8
#> [3,]    3    6    9
  • Apply 함수는 행렬의 행 또는 열 방향으로 특정 함수를 적용한다.
  • apply(array, 방향, 함수)
  • 1: 행, 2: 열
apply(x, 1, function(x) {2*x}) 
#>      [,1] [,2] [,3]
#> [1,]    2    4    6
#> [2,]    8   10   12
#> [3,]   14   16   18
# apply 함수는 vector 에 적용할 수 없다. 
# dim attribute 가 있어야 함
# apply(c(1,2,3), 1, function(x) {2*x}) 
# 에러 출력

lapply

  • apply 함수의 단점은 input 으로 array 만 입력할 수 있다는 것이다. 
  • 일반적으로 vector 를 input 넣는 경우가 많은데, 이를 위해 lapply 가 존재한다. 
  • 입력으로 vector 또는 list 를 받아 list 를 반환한다.
result <- lapply(1:3, function(x) x*2)
result
#> [[1]]
#> [1] 2
#> 
#> [[2]]
#> [1] 4
#> 
#> [[3]]
#> [1] 6
unlist(result) # vector 로 바꾸고 싶으면 unlist 
#> [1] 2 4 6
  • 데이터 프레임에도 lapply 를 적용할 수 있다.
  • 데이터 프레임은 list 에 기반한 s3 object 이기 때문
head(iris)
#>   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#> 1          5.1         3.5          1.4         0.2  setosa
#> 2          4.9         3.0          1.4         0.2  setosa
#> 3          4.7         3.2          1.3         0.2  setosa
#> 4          4.6         3.1          1.5         0.2  setosa
#> 5          5.0         3.6          1.4         0.2  setosa
#> 6          5.4         3.9          1.7         0.4  setosa
typeof(iris)
#> [1] "list"
lapply(iris[, 1:4], mean)
#> $Sepal.Length
#> [1] 5.843333
#> 
#> $Sepal.Width
#> [1] 3.057333
#> 
#> $Petal.Length
#> [1] 3.758
#> 
#> $Petal.Width
#> [1] 1.199333
y <- lapply(iris[, 1:4], function(x) {x > 3})
head(lapply(y, function(x) x[1:5]))
#> $Sepal.Length
#> [1] TRUE TRUE TRUE TRUE TRUE
#> 
#> $Sepal.Width
#> [1]  TRUE FALSE  TRUE  TRUE  TRUE
#> 
#> $Petal.Length
#> [1] FALSE FALSE FALSE FALSE FALSE
#> 
#> $Petal.Width
#> [1] FALSE FALSE FALSE FALSE FALSE

sapply

  • sapply 는 list 대신 행렬 or 벡터로 반환한다.
  • lapply는 list 를 반환하므로 list 를 다시 unlist 하는 것이 번거롭다.
  • lapply 의 wrapper 이다.
y <- sapply(iris[,1:4], function(x) {x > 3})
typeof(y) # Logical matrix 로 반환한다.
#> [1] "logical"
class(y)
#> [1] "matrix"
head(y)
#>      Sepal.Length Sepal.Width Petal.Length Petal.Width
#> [1,]         TRUE        TRUE        FALSE       FALSE
#> [2,]         TRUE       FALSE        FALSE       FALSE
#> [3,]         TRUE        TRUE        FALSE       FALSE
#> [4,]         TRUE        TRUE        FALSE       FALSE
#> [5,]         TRUE        TRUE        FALSE       FALSE
#> [6,]         TRUE        TRUE        FALSE       FALSE
z <- sapply(iris[, 1:4], mean)
z
#> Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
#>     5.843333     3.057333     3.758000     1.199333

vapply

  • FUN.VALUE argument 에 output format 을 명확히 정의해서 더 안전함
y <- vapply(iris[, 1:4], function(x) {x > 3}, numeric(length(iris[, 1]))) # numeric vector 로의 반환 
typeof(y)
#> [1] "double"
class(y)
#> [1] "matrix"
head(y)
#>      Sepal.Length Sepal.Width Petal.Length Petal.Width
#> [1,]            1           1            0           0
#> [2,]            1           0            0           0
#> [3,]            1           1            0           0
#> [4,]            1           1            0           0
#> [5,]            1           1            0           0
#> [6,]            1           1            0           0

z <- vapply(iris[, 1:4], function(x) {x > 3}, logical(length(iris[, 1]))) # logical vector 로의 반환 
typeof(z)
#> [1] "logical"
class(z)
#> [1] "matrix"
head(z)
#>      Sepal.Length Sepal.Width Petal.Length Petal.Width
#> [1,]         TRUE        TRUE        FALSE       FALSE
#> [2,]         TRUE       FALSE        FALSE       FALSE
#> [3,]         TRUE        TRUE        FALSE       FALSE
#> [4,]         TRUE        TRUE        FALSE       FALSE
#> [5,]         TRUE        TRUE        FALSE       FALSE
#> [6,]         TRUE        TRUE        FALSE       FALSE

tapply

  • tapply 는 그룹별 처리를 위한 함수이다. 
  • 그룹을 인자로 주고 (factor 형), 원소별 처리가 아니라 그룹별 처리를 함 
tapply(1:10, rep(c(1,2), each=5), sum) # 그룹 인자는 알아서 factor 형으로 변환이 됨 
#>  1  2 
#> 15 40
  • iris 데이터에서 Species 별로 Sepal.Length 의 평균 구하기
tapply(iris$Sepal.Length, iris$Species, mean)
#>     setosa versicolor  virginica 
#>      5.006      5.936      6.588

mapply

  • sapply 와 비슷하지만 여러개의 인자를 넘긴다.
mapply(function(i, s) {
       sprintf(" %d%s ", i, s)
}, 1:3, c("a", "b", "c"))
#> [1] " 1a " " 2b " " 3c "
  • iris 데이터에서 Sepal.Length 와 Sepal.Width 를 합친 새로운 변수 생성
  • 인자가 몇 개 올지 모르므로, 앞선 함수들에서는 data 를 먼저 인자로 넘겼지만 mapply 에서는 함수를 먼저 인자로 넘긴다.
y <- mapply(`+`, iris$Sepal.Length, iris$Sepal.Width)
# 위 코드는 mapply(function(a, b) a+b, iris$Sepal.Length, iris$Sepal.Width) 와 같다. 
head(iris[,c("Sepal.Length", "Sepal.Width")])
#>   Sepal.Length Sepal.Width
#> 1          5.1         3.5
#> 2          4.9         3.0
#> 3          4.7         3.2
#> 4          4.6         3.1
#> 5          5.0         3.6
#> 6          5.4         3.9
y[1:5]
#> [1] 8.6 7.9 7.9 7.7 8.6


  • 행인 2020.05.12 17:16

    코드 가독성이 너무 떨어지는군요.

    • 행인2 2020.11.27 13:09

      눈이 침침하신가봐요 안경 쓰세요!

    • 행인3 2020.12.31 11:01

      아니면 라섹이라도 ㅎㅎ

    • 행인4 2021.01.05 08:49

      두세줄짜리 코드를 못 읽으면 어떡하자는 건지?

    • 행인5 2021.02.02 11:54

      저것도 읽기 어려우면 일찍이 프로그래밍은 멀리 하세요 ㅋ

  • 행인 2021.05.10 23:18

    돋보기 가져 왔습니다. 코드 가독성이 좋군요

  • YOoOJM 2021.06.17 09:17 신고

    ㅋㅋㅋㅋ웃기네 글씨가 얇아서 그러신듯


다양한 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



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을 이용할 수 있습니다. 

  • 귀곡자 2020.06.02 14:18

    잘 배우고 갑니다

  • R 2020.10.17 16:16

    감사합니다~

  • 아아 2021.04.13 23:49

    감사합니다. 근데 새 버전 설치하고 나면, 이전 버전 폴는 지워도 되나요??

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 분포를 이용한다. 

  • bsi1209 2020.05.18 16:05

    도움 많이 됐습니다. 감사합니다.

R igraph 설치 오류 해결


문제 : libgfortran.so.4 오류로 인해 igraph 패키지가 설치되지 않음


해결 : 

sudo apt-get install gcc-4.8-base libgfortran4


library(devtools)

install_github("igraph/rigraph")

R Default library path 바꾸기 (공용 라이브러리 경로 설정)


.libPaths() 를 하면 현재 사용하고 있는 라이브러리 경로들이 나온다.


[1] "C:/Users/사용자이름/Documents/R/win-library/3.3"

[2] "C:/Program Files/R/R-3.3.2/library"


R은 여기서 첫 번째 element를 Default로 쓰는데, 이것은 보통 개인 폴더로 잡혀 있다. 

하지만 다중 사용자 환경에서는 내가 설치한 라이브러리를 남이 사용하도록 해야할 경우도 있는데 이때는 아래와 같이 하면 된다. 


1. R_LIBS 환경변수를 공용 라이브러리 경로로 설정함 

2. R에서 Sys.getenv("R_LIBS")를 통해 환경변수가 잘 잡혔는지를 볼 수 있다.

3. .libPaths()를 다시 입력하면 공용 라이브러리 경로가 1번으로 잡혀있는 것을 확인할 수 있다. 

 

[1] "C:/Program Files/R/R-3.3.2/library" 

[2] "C:/Users/사용자이름/Documents/R/win-library/3.3"


이후에는 패키지를 설치하면 모든 사용자가 이용할 수 있게 된다. 

  • JHY 2020.09.05 02:21

    질문이 하나 있습니다.

    [1] "C:/Users/사용자이름/Documents/R/win-library/3.3"

    [2] "C:/Program Files/R/R-3.3.2/library"


    R의 버전이 다르긴 하지만, 위 2 경로 중에 한 곳에만 설치되어 있어도 되는 것인가요?
    아니면 두곳 다 설치되어 있어야 정상인 것인가요?

    감사합니다.

  • spogood744 2020.11.24 01:04

    유용한 내용 정말 잘 보고 갑니다


PBC 데이터를 통한 생존분석


PBC 데이터 다운받을 수 있는 사이트

http://www4.stat.ncsu.edu/~boos/var.select/pbc.html

http://www4.stat.ncsu.edu/~boos/var.select/pbc.dat.txt


이 데이터셋은 mayo clinic에서 수집된 데이터로 일차성 담즙성 간경화증(primary biliary cirrhosis) 환자들의 생존에 대한 데이터입니다. 


데이터 핸들링

library(survival)

data <- read.table("pbc.dat")

# Column rename
colnames(data) <- c('id', 'futime', 'status', 'drug', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'edema', 'bili', 'chol', 'albumin', 'copper', 'alk_phos', 'sgot', 'trig', 
                   'platelet', 'protime', 'stage')

head(data)


id
futime
status
drug
age
sex
ascites
hepato
spiders
edema
bili
chol
albumin
copper
alk_phos
sgot
trig
platelet
protime
stage
1400212146411111.014.52612.601561718.0137.9517219012.24
24500012061710110.01.13024.14547394.8113.528822110.63
31012212559400000.51.41763.48210516.096.105515112.04
41925211999410110.51.82442.54646121.860.639218310.34
51504121391810110.03.42793.53143671.0113.157213610.93
62503222420110100.00.82483.9850944.093.0063.11.03

생존분석 자료의 기본 자료 구성


1. Survival 데이터에서는 time을 나타내는 변수와 event를 나타내는 변수가 있다. 

2. time은 event 까지의 시간을 나타내며, 

3. event 변수는 우리의 관심 event (예를 들어, 사망) 혹은 censoring을 나타낸다. (이는 주로 0, 1로 나타내어진다.)


PBC 데이터에서는  status = 0=alive, 1=liver transplant, 2=dead 인데,  0,1을 중도 절단된 censoring, 2 를 관심 event로 코딩을 해봅시다.

data$status[(data$status == 1)] = 0
data$status[(data$status == 2)] = 1


drug = 1 : 페니실린

drug = 2 : 플라시보이므로 , 해석상의 이점을 얻기 위해

drug = 1 페니실린, drug = 0은 플라시보로 재코딩해줍니다.


data <- data[data$drug != '.', ]
data$drug <- as.character(data$drug)
data$drug[(data$drug == 2)] = 0
data$drug <- factor(data$drug)


1. KM estimation


우선 Survival function을 추정하는 비모수적인 방법으로 많이 쓰이는 Kaplan-Meier analysis 를 해보겠습니다.


# event, censoring 구분
Y = Surv(data$futime, data$status)

# KM estimation
fit = survfit(Y~data$drug)
summary(fit)


우선 Surv라는 함수를 통해 무엇이 time과 event인지를 알려줍니다. 그리고 drug에 따라 생존을 분석할 것이므로 ~data$drug를 입력하면 아래처럼 lifetable을 만들어줍니다.


Call: survfit(formula = Y ~ data$drug)

                data$drug=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   51    154       1    0.994 0.00647        0.981        1.000
   77    153       1    0.987 0.00912        0.969        1.000
  110    152       1    0.981 0.01114        0.959        1.000
  130    151       1    0.974 0.01282        0.949        0.999
  186    150       1    0.968 0.01428        0.940        0.996
  191    149       1    0.961 0.01559        0.931        0.992
  207    148       1    0.955 0.01679        0.922        0.988
  216    147       1    0.948 0.01788        0.914        0.984
  264    146       2    0.935 0.01986        0.897        0.975
  304    144       1    0.929 0.02075        0.889        0.970
  321    143       1    0.922 0.02160        0.881        0.965
  326    142       1    0.916 0.02240        0.873        0.961
  460    141       1    0.909 0.02317        0.865        0.956
  549    140       1    0.903 0.02389        0.857        0.951
  552    139       1    0.896 0.02459        0.849        0.946
  597    138       1    0.890 0.02525        0.841        0.941
  611    137       1    0.883 0.02589        0.834        0.935
  708    136       1    0.877 0.02650        0.826        0.930
  733    135       1    0.870 0.02709        0.819        0.925
  769    134       1    0.864 0.02765        0.811        0.920
  786    133       1    0.857 0.02820        0.804        0.914
  790    131       1    0.851 0.02873        0.796        0.909
  797    130       1    0.844 0.02925        0.789        0.903
  850    128       1    0.837 0.02975        0.781        0.898
  853    127       1    0.831 0.03024        0.774        0.892
  859    126       1    0.824 0.03071        0.766        0.887
  890    125       1    0.818 0.03116        0.759        0.881
  930    124       1    0.811 0.03160        0.751        0.875
  943    123       1    0.804 0.03203        0.744        0.870
  974    122       1    0.798 0.03244        0.737        0.864
 1080    118       1    0.791 0.03286        0.729        0.858
 1165    115       1    0.784 0.03328        0.722        0.852
 1212    114       1    0.777 0.03370        0.714        0.846
 1217    111       1    0.770 0.03411        0.706        0.840
 1356    103       1    0.763 0.03459        0.698        0.834
 1413    101       1    0.755 0.03506        0.690        0.827
 1427     98       1    0.748 0.03554        0.681        0.821
 1444     95       1    0.740 0.03603        0.672        0.814
 1487     93       1    0.732 0.03651        0.664        0.807
 1536     91       1    0.724 0.03698        0.655        0.800
 1786     79       1    0.715 0.03763        0.645        0.792
 1847     76       1    0.705 0.03829        0.634        0.784
 2090     69       1    0.695 0.03908        0.622        0.776
 2419     56       1    0.683 0.04030        0.608        0.766
 2466     53       1    0.670 0.04155        0.593        0.756
 2503     51       1    0.657 0.04276        0.578        0.746
 2769     40       1    0.640 0.04473        0.558        0.734
 2796     38       1    0.623 0.04662        0.538        0.722
 2847     35       1    0.605 0.04857        0.517        0.709
 3090     32       1    0.587 0.05060        0.495        0.695
 3170     29       1    0.566 0.05275        0.472        0.680
 3244     28       1    0.546 0.05460        0.449        0.664
 3358     26       1    0.525 0.05640        0.425        0.648
 3395     24       1    0.503 0.05814        0.401        0.631
 3428     22       1    0.480 0.05983        0.376        0.613
 3445     21       1    0.457 0.06119        0.352        0.595
 3762     15       1    0.427 0.06427        0.318        0.573
 3839     13       1    0.394 0.06719        0.282        0.551
 3853     12       1    0.361 0.06916        0.248        0.526

                data$drug=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   41    158       1    0.994 0.00631        0.981        1.000
   71    157       1    0.987 0.00889        0.970        1.000
  131    156       1    0.981 0.01086        0.960        1.000
  140    155       1    0.975 0.01250        0.950        0.999
  179    154       1    0.968 0.01393        0.941        0.996
  198    153       1    0.962 0.01521        0.933        0.992
  223    152       1    0.956 0.01637        0.924        0.988
  334    151       1    0.949 0.01744        0.916        0.984
  348    150       1    0.943 0.01844        0.908        0.980
  388    149       1    0.937 0.01937        0.900        0.975
  400    148       1    0.930 0.02025        0.892        0.971
  515    147       1    0.924 0.02108        0.884        0.966
  673    145       1    0.918 0.02187        0.876        0.962
  694    144       1    0.911 0.02263        0.868        0.957
  750    141       1    0.905 0.02337        0.860        0.952
  762    140       1    0.898 0.02408        0.852        0.947
  799    139       1    0.892 0.02476        0.845        0.942
  824    138       1    0.885 0.02541        0.837        0.937
  904    134       1    0.879 0.02607        0.829        0.931
  971    132       1    0.872 0.02671        0.821        0.926
  980    131       1    0.866 0.02732        0.814        0.921
  999    130       1    0.859 0.02791        0.806        0.915
 1000    129       1    0.852 0.02848        0.798        0.910
 1012    128       1    0.846 0.02902        0.791        0.904
 1037    127       1    0.839 0.02955        0.783        0.899
 1077    126       1    0.832 0.03005        0.775        0.893
 1083    125       1    0.826 0.03054        0.768        0.888
 1152    124       1    0.819 0.03101        0.760        0.882
 1170    122       1    0.812 0.03148        0.753        0.876
 1191    121       2    0.799 0.03236        0.738        0.865
 1235    117       1    0.792 0.03279        0.730        0.859
 1297    114       1    0.785 0.03323        0.723        0.853
 1350    111       1    0.778 0.03368        0.715        0.847
 1360    110       1    0.771 0.03411        0.707        0.841
 1434    105       1    0.764 0.03456        0.699        0.834
 1492    100       1    0.756 0.03505        0.690        0.828
 1576     97       1    0.748 0.03554        0.682        0.821
 1657     93       1    0.740 0.03606        0.673        0.814
 1682     92       1    0.732 0.03655        0.664        0.807
 1690     91       2    0.716 0.03748        0.646        0.793
 1741     87       1    0.708 0.03794        0.637        0.786
 1827     82       1    0.699 0.03845        0.628        0.779
 1925     78       1    0.690 0.03899        0.618        0.771
 2055     72       1    0.681 0.03960        0.607        0.763
 2081     71       1    0.671 0.04019        0.597        0.755
 2105     70       1    0.661 0.04074        0.586        0.746
 2224     65       1    0.651 0.04137        0.575        0.738
 2256     63       1    0.641 0.04198        0.564        0.729
 2288     61       1    0.630 0.04259        0.552        0.720
 2297     60       1    0.620 0.04315        0.541        0.710
 2386     54       1    0.608 0.04385        0.528        0.701
 2400     53       1    0.597 0.04450        0.516        0.691
 2540     47       1    0.584 0.04533        0.502        0.680
 2583     42       1    0.570 0.04634        0.486        0.669
 2598     41       1    0.556 0.04725        0.471        0.657
 2689     38       1    0.542 0.04822        0.455        0.645
 3086     28       1    0.522 0.05023        0.433        0.631
 3222     24       1    0.501 0.05264        0.407        0.615
 3282     22       1    0.478 0.05495        0.381        0.599
 3574     18       1    0.451 0.05795        0.351        0.580
 3584     17       1    0.425 0.06032        0.322        0.561
 4079      8       1    0.372 0.07247        0.254        0.545
 4191      7       1    0.319 0.07922        0.196        0.519


par(mai=c(1,1,1,1))
plot(fit,main="KM Curve", xlab="Time(Week)", ylab=expression(paste(hat(S),"(t)")))

# Log rank test (difference between groups)
log_rank = survdiff(Surv(data$futime, data$status) ~ data$drug, data)
log_rank

결과를 보시면 플라시보와 페니실린 처방군간에 생존률에 큰 차이가 없는 것을 확인할 수 있습니다.


두 군간에 차이가 있는지 검정하는 방법인 log_rank 테스트에서도 p value > 0.05로 두 군간에 차이가 유의하지 않습니다.


Call:
survdiff(formula = Surv(data$futime, data$status) ~ data$drug, 
    data = data)

              N Observed Expected (O-E)^2/E (O-E)^2/V
data$drug=0 154       60     61.8    0.0513     0.102
data$drug=1 158       65     63.2    0.0502     0.102
 

Chisq= 0.1 on 1 degrees of freedom, p= 0.7



2. Cox regression


생존분석을 다변수로 할 때, Cox regression을 주로 활용합니다. 이는 Cox proportional hazard model로 불리기도 합니다.


Cox regression은 다음과 같은 이유로 많이 쓰입니다.


1. 다양한 변수들의 계수와 hazard ratio를 추정하는데 좋다. -> 이를 통해 변수별 생존에 대한 중요도를 비교하기 쉽다.

2. Cox regression은 robust 하다. (참에 매우 근접한 값을 추정해준다.)

3. 계산된 hazard가 양수이기 때문에 타당한 값이다.

4. baseline hazard를 몰라도 회귀계수, hazard ratio를 추정할 수 있다.


Cox regression의 이론적 배경에 대해서는 다음에 기회가 있을 때 다루어보도록 하겠습니다. 


R에서 Cox regression을 수행하는 것은 매우 간단한데, 앞서 정의한 event, time 변수인 Y를 종속변수로 놓고 분석대상인 변수들을 ~ +로 연결하면 됩니다.  

# Cox regression
cox = coxph(Y~data$drug + data$age + data$sex + data$stage + data$edema + data$bili + data$albumin)
summary(cox) 
Call:
coxph(formula = Y ~ data$drug + data$age + data$sex + data$stage + 
    data$edema + data$bili + data$albumin)

  n= 312, number of events= 125 

                   coef  exp(coef)   se(coef)      z Pr(>|z|)    
data$drug1    2.569e-02  1.026e+00  1.855e-01  0.138 0.889882    
data$age      6.736e-05  1.000e+00  2.651e-05  2.541 0.011057 *  
data$sex     -5.940e-01  5.521e-01  2.550e-01 -2.330 0.019832 *  
data$stage1  -2.248e+00  1.056e-01  1.023e+00 -2.197 0.027992 *  
data$stage2  -8.573e-01  4.243e-01  2.950e-01 -2.906 0.003659 ** 
data$stage3  -4.884e-01  6.136e-01  2.186e-01 -2.234 0.025487 *  
data$stage4          NA         NA  0.000e+00     NA       NA    
data$edema    9.740e-01  2.649e+00  3.148e-01  3.095 0.001971 ** 
data$bili     1.263e-01  1.135e+00  1.529e-02  8.263  < 2e-16 ***
data$albumin -9.288e-01  3.950e-01  2.556e-01 -3.634 0.000279 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

             exp(coef) exp(-coef) lower .95 upper .95
data$drug1      1.0260     0.9746   0.71322    1.4760
data$age        1.0001     0.9999   1.00002    1.0001
data$sex        0.5521     1.8112   0.33497    0.9101
data$stage1     0.1056     9.4656   0.01423    0.7844
data$stage2     0.4243     2.3568   0.23801    0.7564
data$stage3     0.6136     1.6296   0.39979    0.9419
data$stage4         NA         NA        NA        NA
data$edema      2.6486     0.3776   1.42921    4.9083
data$bili       1.1347     0.8813   1.10117    1.1692
data$albumin    0.3950     2.5315   0.23936    0.6519

Concordance= 0.839  (se = 0.029 )
Rsquare= 0.427   (max possible= 0.983 )
Likelihood ratio test= 174  on 9 df,   p=<2e-16
Wald test            = 186.2  on 9 df,   p=<2e-16
Score (logrank) test = 287  on 9 df,   p=<2e-16


그러면 이렇게 회귀계수와 회귀 계수의 exponential을 돌려주게 됩니다. Cox regression이 좋은 점은 회귀계수에 exponential을 취했을 때, 해당 변수가 1단위 증가했을 때의 hazard ratio가 됩니다. 예를 들어, age의 경우, 다른 변수들을 보정하고, 1살 증가했을 때, hazard가 1.0001배 증가한다. 라고 볼 수 있습니다. sex의 경우, 다른 모든 변수를 보정했을 때, 남자에 비해 여자가 hazard가 0.55 배라고 볼 수 있습니다 (이 데이터에서 0=male, 1=female로 코딩이 되어있기 때문에). p-value를 보면 어떠한 변수가 생존에 영향을 크게 끼치는 지를 비교해볼 수 있습니다.  



Jupyter notebook 에서 R을 이용하기


데이터 분석을 할 때, PythonR을 동시에 이용하시는 분들이 많습니다.

저 같은 경우도 머신러닝, 딥러닝 등을 할 때는 Python을, 데이터 처리, 가공, 통계분석을 할 때는 R을 선호하는데요. 이 경우 Jupyter notebook 등의 Python 개발 환경과 R studio 를 번갈아가면서 이용해야해서 불편한 점이 많았습니다. 


본 포스팅에서는 Jupyter notebook이라는 하나의 개발 환경에서 R과 Python을 같이 이용할 수 있게 해주는 IRkernel을 설치하는 것을 포스팅해보았습니다.


IRkernel을 이용하면 Jupyter notebook에서도 현재 컴퓨터에 설치된 R 커널을 이용할 수 있게 해줍니다.

IRkernel은 R 커널을 Jupyter notebook에서 이용할 수 있게 해주는 패키지입니다. 

https://github.com/IRkernel/IRkernel


우선 당연히 선행적으로 Jupyter가 설치되어있어야하며,


주의할 점은 conda 가상환경 쓰시는 분들은

이걸 jupyter notebook이 설치된 conda 가상환경을 키고 IRkernel 설치를 진행하여야합니다.


먼저 devtools를 설치합니다.

install.packages('devtools')


다음 devtools를 이용해 IRkernel을 설치합니다.


devtools::install_github('IRkernel/IRkernel')

# or devtools::install_local('IRkernel-master.tar.gz')

IRkernel::installspec()  # to register the kernel in the current R installation



이렇게 R kernel이 생기게 됩니다.



저는 R version 3.5.1을 이용중이었는데 동일한 버전을 Jupyter notebook에서도 이용할 수 있습니다.


---- 추가 (2020-09-28)

conda install r-irkernel

위 conda 커맨드로 anaconda 가상환경 밑에 R 과 IRKernel 을 동시에 설치할 수도 있습니다. 


특정 R 버전에 설치하기 위해서는 아래 참고 (2020-10-16 추가)

http://salvatoregiorgi.com/blog/2018/10/16/installing-an-older-version-of-r-in-a-conda-environment/

  • 데이터사이언티스트꿈나무 2019.07.29 21:52

    안녕하세요. R사용하다가 현재는 파이썬 공부중인 학생입니다. 요즘 파이썬 공부하면서 느낀점은 파이썬의 datetiem객체가 너무 불편하다는 점입니다. R lubridate는 정말 편한 반면에요. 그 밖에도 pandas 문법이 tidyverse보다 불편하기도 하고 등등 머신러닝/딥러닝 제외하면 아무리 생각해도 R이 편한 것 같습니다.

    • Deepplay 2019.08.02 16:19 신고

      네 동의하는 바입니다. 저도 tabular 형식 데이터 프로세싱과 통계 분석을 할 때는 R 을 주로 사용하고 있습니다. 하지만 분석 내용 토대로 어플리케이션을 구축하는 목적에서는 일반 프로그래밍 언어인 파이썬이 강점이 있는것 같습니다.

  • QA 2019.08.05 21:28

    R과 파이썬을 선택하려고하면 어떤것을 선택하고 시작하는것을 추천드리나요? 전산통계전공이고 파이썬과 R 등 프로그래밍경험이 있고 통계과목들을 수강한 베이스입니다.

    • Deepplay 2019.08.06 17:50 신고

      안녕하세요. R과 Python 중에 무엇을 선택하냐는 문제는 이 언어들을 통해 무엇을 하느냐에 따라 다르다고 생각합니다. 논문, 프로젝트 등을 위한 데이터처리/통계 분석에는 R 이 강점이 있고, 딥러닝/머신러닝 모형 구축/어플리케이션 구현에는 Python 이 강점이 있습니다.

      전산통계 전공이시고, 데이터 과학 관련 공부, 프로젝트가 목적이라면 R로 시작하는 것을 추천드립니다. R 이 데이터과학분야의 전통 언어라고 볼 수 있어 기본 함수 (base, stat) 가 잘되어있고, 데이터처리 패키지에 있어서는 강력한 패키지 (tidyverse 등) 들이 많기 때문에 시작하기 쉽습니다. R을 일정수준이상 다룰 수 있으면 Python 도 쉽게 할 수 있다고 생각합니다.

  • python and R 2019.08.13 16:34

    rstudio 안에서 rpython을 이용하는 방식처럼 jupyter notebook 안에서 Rpy2을 쓰면서 서로 같이 쓰시는건가요??! 아니면 R이랑 python 스크립트 창을 따로 만들어야하는데 어떻게 같이 번갈아가면서 쓰는지 궁금합니다.

    • Deepplay 2019.08.17 13:32 신고

      안녕하세요. R과 Python 을 같이 쓴다는 것이 같은 개발환경에서 R과 Python 을 같이 쓴다는 의미였습니다! 제 경우에는 rpy2 는 현재는 쓰지 않고 있고, R 커널을 실행할 수 있는 스크립트 창을 따로만들어서 쓰고 있습니다.