broomパッケージを使ってみる

broomパッケージの上手な使い方を考えてる。

library(stargazer)
library(broom)
library(dplyr)
library(ggplot2)

res1 <- lm(yield ~ block + N + P + K, npk)
res2 <- lm(yield ~ block + N + P * K, npk)

# res <- res1
res <- res2

anova(res)
res %>% 
  aov() %>% 
  tidy()
summary(res)
tidy(res)

# 多重比較
res %>% 
  aov() %>% 
  TukeyHSD() %>% 
  tidy() %>%
  mutate(pair = paste(term, comparison, sep="_")) %>%
  ggplot(aes(colour=cut(adj.p.value, c(0, 0.01, 0.05, 1),
                        label=c("p<0.01","p<0.05","Non-Sig")))) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(pair, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_point(aes(pair, estimate)) +
  labs(colour="significance", title = "Tukey Honest Significant Differences") +
  coord_flip()

glance(res)
augment(res)

res %>% 
  tidy(conf.int = TRUE) %>% 
  filter(term != "(Intercept)") %>% 
  ggplot(aes(estimate, term, color = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 0)

stargazer(res1, res2, type="text", column.labels = c("res1","res2"),) 

続々・複数のファイルを一度に読みこむ方法

過去2回の続きでこれが最後(多分)。
複数のファイルを一度に読みこむ方法 - bob3’s blog
続・複数のファイルを一度に読みこむ方法 - bob3’s blog

前回、data.table::fread() を foreach で並列化したパターンを追加したのだけど、data.table::fread() はそもそも並列化されているということを思い出しました。
大変な無駄をしていたわけです。

このままだと並列化によって速くなるか否かが不明なままですので、readr::read_csv() を並列化して速さを確認してみます。
結論から言うとやっぱり fread & rbindlist が最速。
でも read_csv() も並列化した分は速くなってます。

f:id:bob3:20190709231809p:plain
read_csvの並列化追加

まあ、通常は fread & rbindlist でいいと思います。

# demo data ---------------------------------------------------------------
library(tidyverse)
library(data.table)
library(doParallel)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 365
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}

# read_csv + map_dfr -----------------------------------------------------------------
read_csv_and_map_dfr <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(
      read_csv,
      progress = FALSE,
      .id = "file_name"
    )
}

# fread + rbindlist ---------------------------------------------------------
fread_and_rbindlist <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- map(file_names, fread, showProgress = FALSE)
  setattr(tmp, "names", file_names)
  rbindlist(tmp, idcol = "file_name")
}

# 並列化 ---------------------------------------------------------------------
cores <- detectCores(logical = FALSE) # コア数確認
cl <- makeCluster(cores)
registerDoParallel(cl)
clusterEvalQ(cl, c(library(readr), library(foreach)))

my_foreach <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp1 <- foreach(i = file_names) %dopar% {
    read_csv(i)
  }
  bind_rows(tmp1, .id = "files") %>%
    left_join(tibble(file_names = file_names, files = as.character(1:length(file_names))), by = "files") %>%
    select(-files)
}

# microbenchimark ---------------------------------------------------------

gc()
gc()

compare <- microbenchmark(
  "read_csv & map_dfr" = read_csv_and_map_dfr(),
  "fread & rbindlist" = fread_and_rbindlist(),
  "doParallel & read_csv & bind_rows" = my_foreach(),
  times = 30L
)

compare

autoplot(compare) + labs(title = "R 3.6.1")

続・複数のファイルを一度に読みこむ方法

前々回の続編。

前々回の記事を書いた後の反応にこんなのがありました。

なるほど、並列化という手があったか。

ということで、まだ動作が不安定な vroom を外して、

  • Tidyverse流(read_csv & map_dfr)
  • data.table流(fread & rbindlist)
  • data.table流並列化(fread & rbindlist & foreach)

の3パターンで試します。

デモデータの用意

必要なパッケージを読み込みつつ、実験用のデモデータの用意。

library(tidyverse)
library(data.table)
library(doParallel)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 365 # 前回は不安定なvroomに合わせて10ファイルだったのを増やした。
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}

Tidyverse流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(
      read_csv,
      progress = FALSE,
      .id = "file_name"
    ) -> dat

data.table流

  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- map(file_names, fread, showProgress = FALSE)
  setattr(tmp, "names", file_names)
  dat <- rbindlist(tmp, idcol = "file_name")

data.table流並列化

cores <- detectCores(logical = FALSE) # コア数確認
cl <- makeCluster(cores)
registerDoParallel(cl)
clusterEvalQ(cl, c(library(data.table), library(foreach)))
file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
tmp <- foreach(i = file_names) %dopar% {
    fread(i, showProgress = FALSE)
  }
setattr(tmp, "names", file_names)
rbindlist(tmp, idcol = "files")

比較

ファイル数が増えると Tidyverse 流では遅く感じてしまう。
やはり、data.table流は速い。
期待した並列化だが、バラツキが大きくて一概に速いとは言えない状況。コア数多めのCPU積んでる場合はいいのかもしれない。
書きやすさのTidyverse、速さのdata.tableというところは変わらない模様。

f:id:bob3:20190709010624p:plain
複数ファイルの読み込み速度比較(並列化入り)

# demo data ---------------------------------------------------------------
library(tidyverse)
library(data.table)
library(doParallel)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 365
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}


# read_csv + map_dfr -----------------------------------------------------------------
read_csv_and_map_dfr <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(
      read_csv,
      progress = FALSE,
      .id = "file_name"
    )
}

# fread + rbindlist ---------------------------------------------------------
fread_and_rbindlist <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- map(file_names, fread, showProgress = FALSE)
  setattr(tmp, "names", file_names)
  rbindlist(tmp, idcol = "file_name")
}

# 並列化 ---------------------------------------------------------------------
cores <- detectCores(logical = FALSE) # コア数確認
cl <- makeCluster(cores)
registerDoParallel(cl)
clusterEvalQ(cl, c(library(data.table), library(foreach)))

my_foreach <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- foreach(i = file_names) %dopar% {
    fread(i, showProgress = FALSE)
  }
  setattr(tmp, "names", file_names)
  rbindlist(tmp, idcol = "files")
}
# stopCluster(cl)

# microbenchimark ---------------------------------------------------------

gc()
gc()

compare <- microbenchmark(
  "read_csv & map_dfr" = read_csv_and_map_dfr(),
  "fread & rbindlist" = fread_and_rbindlist(),
  "doParallel & fread & rbindlist" = my_foreach(),
  times = 10L
)

compare

autoplot(compare) + labs(title = "R 3.6.1")

base と Tidyverse と data.table と

恋しさと せつなさと 心強さと」といえば『ストリートファイターII MOVIE』の挿入歌ですね。

さて、初心者にRを教えるとき、base準拠にすべきか、Tidyverse準拠にすべきか、という議論があります。
これは私も悩むところでありまして、さらに言えば data.table ってのもあります。

んで、まずはこの三つの流派(?)が存在してそれぞれこんな書き方をするんだよ、というのを示せるといいのかなーと思うわけです。

そこで試しにこんな処理を base, Tidyverse, data.tableのそれぞれで書いてみました。

  1. データセット flights の変数をdep_delay、arr_delay、carrier、air_time、distanceに絞る。
  2. いずれかの変数にNAを含む行を削除。
  3. 時速をkm/hで算出。air_timeは分単位なので60で除して時間単位に、distanceはマイル単位なので1.60934を乗じてkm単位に。
  4. carrier毎にdep_delay、arr_delay、Hourly speedの平均と件数を算出。
  5. 件数の降順に並べる。

こんな風に書いてみましたが、それぞれの書き方の違いの雰囲気ぐらいは伝わるかな?

base

library(nycflights13)
flights.df <- na.omit(as.data.frame(flights)[, c("dep_delay", "arr_delay", "carrier", "air_time", "distance")])
flights.df$HourlySpeed <- (flights.df$distance * 1.60934) / (flights.df$air_time / 60)
flights.df2 <- aggregate(
  flights.df[, c("dep_delay", "arr_delay", "HourlySpeed") ],
  list(carrier = flights.df$carrier),
  mean
)
flights.df2$n <- table(flights.df$carrier)
flights.df2[order(flights.df2$n, decreasing = TRUE), ]

Tidyverse

library(nycflights13)
library(tidyverse)
flights %>%
  select(dep_delay, arr_delay, carrier, air_time, distance) %>%
  drop_na() %>%
  mutate(HourlySpeed = (distance * 1.60934) / (air_time / 60)) %>%
  group_by(carrier) %>%
  summarise(
    mean_dep_delay = mean(dep_delay),
    mean_arr_delay = mean(arr_delay),
    mean_HourlySpeed = mean(HourlySpeed),
    n = n()
  ) %>% 
  arrange(desc(n))

data.table

library(nycflights13)
library(data.table)
flights.dt <- na.omit(as.data.table(flights)[, c("dep_delay", "arr_delay", "carrier", "air_time", "distance")])
setorder(flights.dt[
  , HourlySpeed := (distance * 1.60934) / (air_time / 60)
  ][
    j = list(
      mean_dep_delay = mean(dep_delay),
      mean_arr_delay = mean(arr_delay),
      mean_HourlySpeed = mean(HourlySpeed),
      n = .N
      ),
    by = carrier
    ],
  -n) -> flights.dt2
flights.dt2

実行速度は?

data.table圧勝!……かと思いきや、これくらいのデータ量だとTidyverseと大差ないですね。

f:id:bob3:20190707223349p:plain
実行速度比較

複数のファイルを一度に読みこむ方法

Rで複数のファイルを一度に読み込みたい場合がたまにある。
例えば、売り上げのデータが日付ごとに分かれたファイルに収められているような場合。
しょっちゅう発生するわけでもなく、いつも忘れてしまって調べるのに時間がかかるのでここにメモしておく。

ここではそうしたファイルが一つのフォルダに収められているものとする。

方法は大きく三つある

  • Tidyverse流(read_csv & map_dfr)
  • data.table流(fread & rbindlist)
  • vroom流(vroom)

デモデータの用意

必要なパッケージを読み込みつつ、実験用のデモデータの用意。

library(tidyverse)
library(data.table)
library(vroom)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 10
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}

Tidyverse流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(read_csv,
      .id = "file_name"
    ) -> dat

data.table流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map(fread) %>%
    rbindlist(idcol = "file_name") -> dat

vroom流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    vroom(id = "file_name") -> dat

比較

Tidyverseに慣れてる人なら read_csv と map_dfr を組み合わせるのが分かりやすいだろう。
data.table流はとにかく速い。爆速。
vroomはファイル数が多くなるとエラーになる。今後に期待。

f:id:bob3:20190707154441p:plain
複数ファイルの読み込み速度比較

# read_csv + map_dfr -----------------------------------------------------------------
read_csv_and_map_dfr <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(read_csv,
      .id = "file_name"
    )
}

# fread + rbindlist ---------------------------------------------------------
fread_and_rbindlist <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map(fread) %>%
    rbindlist(idcol = "file_name")
}

# vroom -------------------------------------------------------------------
# ファイル数が多いとエラーになる場合がある。
my_vroom <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    vroom(id = "file_name")
}

# microbenchimark ---------------------------------------------------------

gc()
gc()

compare <- microbenchmark(
  "read_csv & map_dfr" = read_csv_and_map_dfr(),
  "fread & rbindlist" = fread_and_rbindlist(),
  "vroom" = my_vroom(),
  times = 30L
)

compare

autoplot(compare) +
  labs(title = "R 3.6.1") 

以上。

R本体のアップデート時に入れるパッケージ

まず最初に必ず入れるパッケージ。
あとは必要になったときに入れれば大丈夫。

install.packages(
  c(
    "tidyverse",
    "data.table",
    "psych",
    "vegan",
    "plotly",
    "devtools",
    "colourpicker", # RStudio Addin
    "ggThemeAssist", # RStudio Addin
    "styler" # RStudio Addin
  ),
  Ncpus = getOption("Ncpus", 4L)
)

Rを使ってExcelファイルを読み込んでcsvファイルを出力するまで

那須野さんの記事が興味深かったので、Rでもやってみた。
takuminasuno.com

library(tidyverse) #データ操作用
library(readxl) #Excelファイルの読み込み用
library(stringi) #文字列操作用

# 列名の指定。空列は"blank"とする。
nameList <- c("age", "blank1", "blank2",
              "total_201905", "male_201905", "female_201905", "blank3",
              "total_201812", "male_201812", "female_201812",
              "total-ja_201812", "male-ja_201812", "female-ja_201812")

# Excelファイルの読み込み。読み込む個所はrangeで指定。
# index(age)の補正。空白の除去と英数の半角化。「歳」の追加。
# 数値の桁数の補正(単位を揃える)。2019年は1万倍、2018年は千倍。
dat <- read_excel("05k2-1.xls", range = "A14:M34", col_names = nameList) %>%
  select(-starts_with("blank")) %>%
  mutate(age = age %>%
           str_replace_all("[[:space:]]", "") %>%
           stri_trans_nfkc() %>%
           if_else(str_detect(., pattern="歳"), ., str_c(., "歳"))
         ) %>%
  mutate_at(vars(contains("2019")), ~.*10000) %>%
  mutate_at(vars(contains("2018")), ~.* 1000)

# 横持ちでcsv出力
write_excel_csv(dat, "output_pivotted.csv")

# 縦持ちのデータに変換
dat2 <- dat %>%
  gather(key="sex_yyyymm", value = "value", -age) %>%
  separate(sex_yyyymm, into = c("sex", "yyyymm"), sep = "_")

# 縦持ちでcsv出力
write_excel_csv(dat2, "output_unpivotted.csv")


まあ、RでもPythonでも慣れてる方で書けばいいのかな。