• Home

R에서 Pipe 연산으로 분석하기

magrittr 패키지가 최근에 많은 화제를 불러 일으키고 있다. 사실 이 패키지는 패키지 자체로 유명세를 탔다고 하기 보다는 다른 유명 패키지가 이 패키지를 사용하게 됨으로써 유명세를 탓고 필자도 현재 이 패키지 때문에 dplyr과 같은 패키지를 자연스럽게 사용하게 되었다.

이 패키지는 유닉스에 있는 파이프(|, >) 연산자와 같은 기능을 아래와 같이 R에서 쓸 수 있게 해주는 연산자이다.

library(magrittr)
library(dplyr)
library(lubridate)

#head(iris)
iris %>% head  
##   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
#iris 데이터에 대해서 group by연산을 수행하고 이 결과를 기반으로 Sepal.Length의 평균과 Petal.Length의 평균을 구한다. 
iris %>% group_by(Species) %>% 
  summarise(mofSL=mean(Sepal.Length), mofPL=mean(Petal.Length)) 
## Source: local data frame [3 x 3]
## 
##      Species mofSL mofPL
## 1     setosa 5.006 1.462
## 2 versicolor 5.936 4.260
## 3  virginica 6.588 5.552
#iris 데이터에서 Species 필드를 제외한 필드들 각각의 max, min 값을 구한다. 
iris %>% select(-Species) %>% 
  apply(2, lambda( x ~ c(max(x), min(x))))
##      Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,]          7.9         4.4          6.9         2.5
## [2,]          4.3         2.0          1.0         0.1
#특정 기간의 날짜에서 일요일만 제거한 날짜의 첫 부분을 출력한다. 
seq(ymd("2014-01-01"), ymd('2014-08-01'), by="1 day")  %>% 
  Filter(lambda(x ~  (x %>% weekdays) != "Sunday"),.) %>% head %>% print
## [1] "2014-01-01 UTC" "2014-01-02 UTC" "2014-01-03 UTC" "2014-01-04 UTC"
## [5] "2014-01-06 UTC" "2014-01-07 UTC"

그런데 단순한 Syntactic Sugar의 역할을 할것만 같았던 이 패키지가 필자의 마음을 흔들어 놓았는지 설명할 필요가 있을 것이다.

위 방식의 장점은 아래와 같다.

  • 코드의 복잡도가 늘어가면 유지 보수나 추후 코드 리딩에 문제가 생길 수 있어 적절한 중간 변수를 만들어두는데 위와 같은 방식은 코드의 복잡도가 리딩에 끼치는 악영향이 기존보다 획기적으로 줄어들기 때문에 빈번한 중간 변수 할당에 따른 메모리 소모를 줄일 수 있다.
  • 코드 작성 방향과 데이터 흐름의 방향 및 생각의 방향이 일치화 되어 코드를 이해하기 쉽다.
  • 결국 dplyr을 쓰게되는데, 이덕분에 data.table문법을 따로 배울 필요가 없어지는데 이는 dplyr에서 알아서 처리해 주기 때문이다.

단점은 아래와 같다.

  • 내가 아닌 이 방식을 이해하지 못하는 다른 사람들이 이 코드를 보게 된다면 당황해 할 수 있다.
  • 위 방식만 고수하기 보다는 기존 방식과 위 방식을 겸용하게 됨으로써 자칫 스파게티 코드가 될 위험이 있다.
  • 결국 dplyr을 쓰게 된다. 나는 data.table 사용자인데 말이다(물론 dplyr는 잘 지원해준다).

앞으로 magrittr 패키지는 폭넓게 사용될 것으로 예상이 되는데 이는 인기 패키지인 dplyr과 앞으로 많은 각광을 받게될 ggvis가 이 방식을 사용하면서 많은 사용자 경험을 이끌 것이기 때문이다. 장점이 있는 패키지임에는 분명하기 때문에 사용자는 점점 많아질 것이다.

내가 쓰는 맥용 .Rprofile

금일 서울대학교 강좌를 진행하면서 필자 맥북의 .Rprofile공유를 요청하셔서 이렇게 올려본다.
아래 스크립트의 앞부분은 링크에서 설명을 하지만 맥에서 기본 그래픽 디바이스인 quartz디바이스의 한글 설정을 해주는 스크립트 이다.

그리고 마지막 라인은 자바 설정으로 KoNLP를 맥에서 사용할때 권장하는 설정 내용이다. 물론 다른 OS에서도 설정해주면 좋다. -Xmx8g부분은 본인의 시스템 메모리를 고려해서 적당하게 설정하면 된다.

stringsAsFactors 부분은 데이터를 읽어들일때 R이 임의로 문자열을 factor로 변환시켜 골탕을 먹이는 일이 종종일어나는데 R임의대로 변환을 시키지 말라는 명령어다.

#https://stat.ethz.ch/pipermail/r-sig-mac/2009-October/006731.html
setHook(packageEvent("grDevices", "onLoad"),
        function(...){
        if(capabilities("aqua"))
            grDevices::quartzFonts(
              sans =grDevices::quartzFont(rep("AppleGothic",4)),
              serif=grDevices::quartzFont(rep("AppleMyungjo",4)))
        grDevices::pdf.options(family="Korea1")
        grDevices::ps.options(family="Korea1")
        }
)
attach(NULL, name = "KoreaEnv")
assign("familyset_hook",
       function() {
            macfontdevs=c("quartz","quartz_off_screen")
            devname=strsplit(names(dev.cur()),":")[[1L]][1]
            if (capabilities("aqua") &&
                devname %in% macfontdevs)
                    par(family="sans")
       },
       pos="KoreaEnv")
setHook("plot.new", get("familyset_hook", pos="KoreaEnv"))
setHook("persp", get("familyset_hook", pos="KoreaEnv"))




options(prompt = "R > ")

options(stringsAsFactors=FALSE)

options(java.parameters=c("-Xmx8g", "-Dfile.encoding=UTF-8"))

Resampling으로 회귀계수 검정

re-sampling 방법으로 회귀모형 파라메터를 검정해보기로 한다. 개인적으로는 리샘플링 방법이 직관적이고 설명이 편한 장점이 있어서 자주 쓰고자 하는 소망을 가지고 있고 금번 포스팅도 그 일환이다. 무엇보다 우리 주변의 컴퓨팅 파워는 놀고 있으니 요즘들어 안쓸 이유가 없는 방법이라 생각한다. …

suppressPackageStartupMessages(library(UsingR))
suppressPackageStartupMessages(library(data.table))
suppressPackageStartupMessages(library(extrafont))

data(father.son)

coeffs <- data.table()

#2만번 re-sampling 
for(i in 1:20000){
  rfheigh <- sample(father.son$fheight)
  rsheigh <- sample(father.son$sheight)

  mdl <- lm(rfheigh ~ rsheigh)
  coeffs <- rbind(coeffs, data.table(inte=mdl$coefficients[1], b0=mdl$coefficients[2]))
}

plot(density(coeffs$b0), main="density plot of b0(H0)")

plot of chunk unnamed-chunk-1

xval <- seq(range(father.son$fheight)[1],range(father.son$fheight)[2], length=100)

rows <- as.data.frame(coeffs[sample(1:nrow(coeffs), 30)] )
plot(father.son, main="H0")

for(i in 1:30){
  lines(xval, rows[i, 'b0'] * xval + rows[i, 'inte'], col = "gray")
}

plot of chunk unnamed-chunk-1

#귀무 가설 공간에서 0.05% 미만으로 나올 극단적인 b0값은?
coeffsdec <- coeffs$b0[order(coeffs$b0, decreasing = T)]
coeffsdec[which.max(cumsum(abs(coeffsdec))/sum(abs(coeffsdec)) > 0.0005)]
## [1] 0.1064
rmdl <- lm(fheight ~ sheight, father.son)

rmdl$coefficients
## (Intercept)     sheight 
##     34.1075      0.4889

이미 b0는 0.4889로 귀무가설 공간에서의 극단적인 값(10000번중에 5번 미만)을 상회하는 값이 나와 b0추정치가 유의미하다는 것을 알 수 있다.