Notes of a Data Scientist

Mémorandum d'un spécialiste des données

質問に対して簡潔に回答する

Amazonでは、社内でプレゼンテーションを行うときにパワーポイントの使用が禁止されていることは有名な話である。 パワーポイントの代わりにメモを作成するそうだが、優れたメモを作成する方法がこちらの tweet にまとまっている。

特に質問を受けたときの回答方法が、シンプルながら参考にしたいと思ったので以下に書く。

If you get a question, reply with one of the four Amazon answers:

1. yes.
2. no.
3. A number.
4. I don't know (and will follow up when I do)

私は質問の回答が冗長になりやすいため、これを徹底して実行してみる。

RでERGMを実装してみた

ネットワーク分析に明るい方であれば、2019年に書かれたこの野心的な記事を覚えている人も多いのではないだろうか。

meana0.hatenablog.com

statnet チームが提供している R のパッケージに頼らずに、Exponential Random Graph Model (ERGM) をできるだけスクラッチで実装しようという壮大な試みである。 実装した関数による推定値は、残念ながら {ergm} で推定した値と乖離してしまったという結びになっている。 その原因解明は多くの読者にとって待ち望んでいたことに違いない。

そこで時間を持て余した筆者が、ゴールデンウィークを使って ERGM を (MPLE を除いて) 一から、かつ {Rcpp} を用いて高速に推定できるよう実装した。*1 以下の repository にコードをパッケージとして公開しているので、興味のある方は遊んでみていただけると幸いである。

github.com

flomariagge ~ edges + triangle という単純なモデルを {myergm}{ergm} で 100 回ずつ推定し、その推定値の分布をプロットしたものが以下の図である。 両者は概ね一致しているように思われる。

先の記事ではなぜ推定がうまくいかなかったのか?

気になるのが、一番最初にあげた記事でなぜ ERGM の推定値が {ergm} のそれと比較して大きく異なってしまったのかである。 その原因が、以下の関数にあることを特定した。

#対数尤度比を計算する関数
loglikelihood_ratio <- function(theta,fit,n){
  sample_length <- length(fit$edges)
  sample_index <- (sample_length-n):sample_length
  sampled_edge_count <- fit$edges[sample_index]
  sampled_triangle_count <- fit$triangle[sample_index]
  theta0 <- c(fit$theta0_edges,fit$theta0_triangle)
  g_obs <- c(count_edges(fit$obs_graph),count_triangle(fit$obs_graph))
  first = (theta - theta0) %*% g_obs
  exp_vec <- c()
  for(i in 1:n){
    sam <- exp((theta - theta0) %*% c(sampled_edge_count[i],sampled_triangle_count[i]))
    exp_vec <- c(exp_vec,sam)
  }
  second <- log(sum(exp_vec)/length(exp_vec))
  loglikelihood_ratio <- first - second
  return(loglikelihood_ratio)
}

sample_index <- (sample_length-n):sample_length とあるので、log likelihood ratio を計算するときに、MCMC で生成されたネットワークの数を M としたとき、 (M - n) 番目から M 番目までの連続した全てのサンプルを使っていることになる。

一方、{ergm} には control.ergm() という関数があり、MCMC をチューニングできるが、そのデフォルトは以下のようになっている。

MCMC.interval = 1024,
MCMC.burnin = MCMC.interval * 16,
MCMC.samplesize = 1024,

まず、1024 * 16 を burnin とする。 その後、MCMC で生成されたグラフ同士の間隔が 1024 となるように (MCMC.interval = 1024)、1024 個のグラフをサンプリングする (MCMC.samplesize = 1024)。 そのサンプリングされた 1024 個のグラフを推定値を得るために用いるのである。

これを踏まえると、最初の記事で推定がうまくいかなかった原因は、MCMC でサンプリングされたグラフ同士の間隔が 0 であり、グラフ同士の serial correlation が強いことに原因がありそうである。

そこで、{myergm} でグラフ同士の相関が強く残るようにサンプリングして推定した場合の推定値の分布を見てみよう。 推定がうまくいかなかったケースを、仮に arabou_ergm と呼ぶことにすると、それは以下のように実装できる。

# Number of simulations
R <- 100

df <- 
  foreach(i = 1:R, .combine = "rbind") %do% {
    # Estimate the parameters.
    theta_arabou_ergm <- myergm::myergm_MCMLE(model = flomarriage ~ edges + triangle, 
                                              seed = i, 
                                              p_one_node_swap = 0, 
                                              p_large_step = 0,
                                              p_invert = 0,
                                              MCMC_samplesize = 1024,
                                              MCMC_interval = 1,
                                              second_round = FALSE)
                              
    # Return the output.
    return(data.frame("iter" = i,"param" = c("edges", "triangle"), "arabou_ergm" = theta_arabou_ergm))
  }

{ergm} で推定したケースと比較したとき、それぞれの分布は以下の通りになる。

arabou_ergm の方が推定値にかなりばらつきが出ていることがわかる。

私が実装した {myergm} は、厳密には {ergm} の実装とは異なる点がある。 例えば、{ergm}MCMC のデフォルトの sampler は "tie-no-tie" といって、最初にリンクが存在する dyad か、リンクが存在しない dyad どちらの状態を変化させるかをそれぞれ確率 1/2 で決定し、リンクを加える/消す。 一方 {myergm} では、tie-no-tie sampler は実装しておらず、代わりに Mele (2017) が提案したような large step を許容した sampler、つまりある確率で 1 回のステップで複数の dyad の状態を変化させるような実装になっている。

ただし、これらの MCMC sampling の手法は推定値に大きく影響を与えるものではなく、大事なのは推定に用いるグラフの sampling の方であったようだ。 他にも ERGM には複雑怪奇なことが多いので*2、それらは今後の検討課題としたい。

References

  • Hummel, R. M., Hunter, D. R., & Handcock, M. S. (2012). Improving simulation-based algorithms for fitting ERGMs. Journal of Computational and Graphical Statistics, 21(4), 920-939.
  • Hunter, D. R., Goodreau, S. M., & Handcock, M. S. (2013). ergm. userterms: A Template Package for Extending statnet. Journal of statistical software, 52(2), i02.
  • Mele, A. (2017). A structural model of dense network formation. Econometrica, 85(3), 825-850.

*1:とはいえ、{ergm} の MCMC sampling の速度には遠く及ばない。どうやらグラフの情報を binary tree structure の形で保存していることが、高速で link を加える/消すことを可能にしているようだ (Hunter, Goodreau, and Handcock, 2013)。

*2:Hummel, Hunter, and Handcock (2012) に従い、{ergm} でも実装されている "partial stepping" なるものを {myergm} にも実装しているが、元論文を読んでもその理論的なモチベーションはよくわからなかった。

Rで疎行列に親しむ

仕事でネットワーク分析に苦しむ機会が増えてきた。 その中で疎行列 (sparse matrix) を取り扱う重要性が高まってきている。 本稿の目的は疎行列の意義と使い方を学び、仲良くなることである。

疎行列とは

疎行列とは、要素のほとんどが 0 であるような行列のことである。

巨大な行列を用いた計算が必要になることがある。 行列が巨大になるときに問題となるのは、行列がメモリを圧迫することである。

例えば、1万×1万の単位行列を、以下のように適当に作成する。

I <- diag(nrow = 10000, ncol = 10000)
pryr::object_size(I)
800 MB

単純な単位行列ではあるが、800 MB ものメモリが割り当てられている。

これが、100万×100万の単位行列を普通に作成すると、以下の通りになる。

I <- diag(nrow = 1000000, ncol = 1000000)
Error: cannot allocate vector of size 7450.6 Gb

たった1つの単位行列を作成するために、7 TB ものメモリーが必要となってしまう。 この行列のほとんどは 0 であるが、それらがわざわざ全て保存される結果、このように巨大な容量を食ってしまう。

ということは、行列の要素がほとんど 0 であるという特性を活かせば、無理なく作成することができるはずである。

そこで疎行列の出番である。RではMatrixパッケージを用いて以下のように作成できる。

I_sp <- Matrix::sparseMatrix(i = 1:1000000, j = 1:1000000)
I_sp[1:3, 1:3]
3 x 3 sparse Matrix of class "ngCMatrix"
          
[1,] | . .
[2,] . | .
[3,] . . |
pryr::object_size(I_sp)
8 MB

100万×100万の単位行列も、わずか8MBで作成できる。 対角要素は全てbooleanで表現されているので、そのおかげでさらにメモリ容量が節約できている。

I_sp[1, 1]
[1] TRUE

これを全て1にすると以下の通りとなり、メモリへの割当が多少増加する。

I_one <- as(I_sp, "dgCMatrix")
I_one[1, 1]
[1] 1
pryr::object_size(I_one)
16 MB

行列演算も通常と同じようにできる。 面白くない例だが、疎行列の単位行列を2つ掛け合わせることは以下のようにできる。

I_sp %*% I_sp

ネットワーク分析における疎行列

グラフで誰と誰がつながっているかを表現するために使われる行列を、隣接行列 (adjacency matrix) と呼ぶ。

例えばノードの数が3つで、1 -- 2 -- 3 のようなつながりのある無向グラフは

と表現できる。つまりGの(i, j)要素は、ノードiとjの繋がりを表現する。

グラフが小さいときはまだ上のように表現できる。 しかしグラフが巨大になるとメモリを圧迫してしまうのは、上で見た例の通りである。 そこで、多くのグラフは sparse であるという特性を活かし、グラフを疎行列で表現することがよく使われる。

Rでネットワーク分析を行うときによく使われるパッケージ群の1つとしてstatnetがある。 例えば、指数ランダムグラフモデル (exponential random graph model; ERGM) を推定するときは、 networkでネットワークオブジェクトを作成して、それに基づき推定を行う。 以下はその一例である。

library(statnet)
data(florentine)
est <- ergm(formula = flomarriage ~ edges + absdiff("wealth") + triangle)

しかしそのネットワークに色々と操作を加えるときはそのネットワークを一度隣接行列に置き換えて処理を行う必要が出てくる。 その際にネットワークオブジェクトを隣接行列に置き換えるのは、例えば以下のようなコードがある。

adj <- network::as.matrix.network.adjacency(flomarriage)
head(adj[, 1:3])
           Acciaiuoli Albizzi Barbadori
Acciaiuoli          0       0         0
Albizzi             0       0         0
Barbadori           0       0         0
Bischeri            0       0         0
Castellani          0       0         1
Ginori              0       1         0

しかしどういうわけか、network::as.matrix.network.adjacencyは密行列 (dense matrix) して隣接行列を返す。 これはネットワークが大きい場合、メモリを大きく圧迫してしまう。

それを回避するためには、まずネットワークをedgelist、つまりつながっているノード同士だけを取り出す。

edgelist <- network::as.edgelist(flomarriage)

そこからこのネットワークを疎行列として表現する。

adj_sp <- Matrix::sparseMatrix(i = edgelist[, 1], j = edgelist[, 2], symmetric = TRUE)
pryr::object_size(adj)
3.58 kB
pryr::object_size(adj_sp)
1.82 kB

小さいネットワークだと対した差ではないが、これがノード数が増えると両者の差は爆発的に広がっていく。

Rcppで疎行列をたしなむ

隣接行列を用いた重たい処理が必要なときは、Rの疎行列をRcppに渡して高速化することができる。 ここではRcppArmadilloを用いてやってみる。 やったことはないが、RcppEigenでも同じようにできるらしい。

その処理の例として、ノードiとjのつながりのうち、iとjが同じ特徴量を持つような隣接行列を作ることを考える。 例えば特徴量として性別を考えたとき、同じ性別同士で実際につながっているリンクのみ含んだ隣接行列を作成したい。

まずは適当にグラフを作成し、疎行列で表現する。

# This seed comes from the result of the 2020 Hawks-Giants Japan Series.
set.seed(264)

# Number of nodes
N <- 1000

# Create an edgelist. I consider an undirected network.
edgelist <- 
  tibble::tibble(tail = 1:N, head = 1:N) %>%
  tidyr::expand(tail, head) %>%
  dplyr::filter(tail < head) %>%
  dplyr::mutate(connected = as.integer(unlist(purrr::rbernoulli(n = nrow(.), p = 0.01)))) %>%
  dplyr::filter(connected == 1)

# Create a sparse adjacency matrix
net <- Matrix::sparseMatrix(i = edgelist$tail, j = edgelist$head, 
                            dims = c(N, N), symmetric = TRUE)
net <- as(object = net, Class = "dsCMatrix")

次にノードの特徴量を適当に作成する。

# Vertex feature like gender.
x <- unlist(purrr::rbernoulli(n = N))

そして今やりたい処理をC++のコードで書く。 RcppArmadilloで疎行列を扱うときはarma::sp_matを用いる。

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::sp_mat get_connected_same_feature_dyads_cpp(const arma::sp_mat& adjmat,
                                                  const arma::vec& feature) {
  int n_nodes = adjmat.n_rows;
  arma::sp_mat output(n_nodes, n_nodes);

  for (int i = 0; i < n_nodes; i++) {
    for (int j = 0; j < n_nodes; j++) {
      if (i != j) {
        if (adjmat(i, j) == 1 & feature[i] == feature[j]) {
          output(i, j) = 1;
        }
      }
    }
  }

  return output;
}

実際に実行してみる。

Rcpp::sourceCpp("src/get_connected_same_feature_dyads_cpp.cpp")
output <- get_connected_same_feature_dyads_cpp(net, x)
net[1, 996] == 1 & x[1] == x[996]
[1] TRUE
output[1, 996]
[1] 1

正しく計算できているように思う。

別にやる必要はないが、同じ処理をRで実装すると以下の通りとなる。

get_connected_same_feature_dyads_R <- function(adjmat, x) {
  n_nodes <- length(x)
  output <- Matrix::sparseMatrix(i = {}, j = {}, dims = c(length(x), length(x)))

  for (i in 1:n_nodes) {
    for (j in 1:n_nodes) {
      if (i != j) {
        if (adjmat[i, j] == 1 & x[i] == x[j]) {
          output[i, j] <- TRUE
        }
      }
    }
  }
  return(output)
}

ベンチマークを行うとRcppの方が圧倒的に速い。Rforを書いた時点で負け確定であるのは内緒である。

time <- rbenchmark::benchmark(
  "R" = {
    output <- get_connected_same_feature_dyads_R(net, x)
  }, 
  "cpp" = {
    output <- get_connected_same_feature_dyads_cpp(net, x)
  },
  replications = 1
)

time
  test replications elapsed relative user.self sys.self user.child sys.child
2  cpp            1   0.016     1.00     0.016    0.000          0         0
1    R            1 136.508  8531.75   136.516    0.004          0         0

Rを使ってより高速に処理するには一工夫必要なことも多いが、変に凝ったコードを書くぐらいならRcppで殴ってしまった方が良い。

終わりに

ぜひみなさんも疎行列に親しみ、ネットワーク分析に苦しみましょう。

broomを使って回帰分析の結果を出力する

いわずもがなであるが、『効果検証入門』は大変勉強になる。 因果推論は大学・大学院のときに学んでいるが、その復習に使える。 また、Rを使ってどのように因果推論を行うかについても有益な情報を与える。

その中で、回帰分析の結果の出力方法について私が知らないやり方があったので、ここで書いてみる。 それがbroomというパッケージを使用した方法である。

通常Rで線形回帰分析を行うときはlm関数を用いる。 例のごとく、irisを用いて適当な線形回帰モデルを推定し、その結果をsummaryで出力する。

reg <- 
  lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data = iris)
summary(reg)
Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, 
    data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.82816 -0.21989  0.01875  0.19709  0.84570 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.85600    0.25078   7.401 9.85e-12 ***
Sepal.Width   0.65084    0.06665   9.765  < 2e-16 ***
Petal.Length  0.70913    0.05672  12.502  < 2e-16 ***
Petal.Width  -0.55648    0.12755  -4.363 2.41e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3145 on 146 degrees of freedom
Multiple R-squared:  0.8586,    Adjusted R-squared:  0.8557 
F-statistic: 295.5 on 3 and 146 DF,  p-value: < 2.2e-16

この方法では、下に余計な情報が出てくるなど見づらい、kableExtraなどで出力できないなど不便なことが多い。

そこでbroom::tidyの出番である。 推定された係数の情報のみを抽出し、tibbleとして出力してくれる。

reg <- 
  lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data = iris) %>%
  broom::tidy()
# A tibble: 4 x 5
  term         estimate std.error statistic  p.value
  <chr>           <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)     1.86     0.251       7.40 9.85e-12
2 Sepal.Width     0.651    0.0666      9.77 1.20e-17
3 Petal.Length    0.709    0.0567     12.5  7.66e-25
4 Petal.Width    -0.556    0.128      -4.36 2.41e- 5

大概の場合、関心のあるのは係数の推定値のみであるから、これでほとんど事足りる。 加えてtibbleとして出力してくれるので、例えばkableExtra::kable(reg)として出力することができる。

論文への出力のためのtableはstargazerを用いることが多いが、書き方によってはこのbroomを用いることで置き換えることが可能かもしれない。

purrrを使ってみる

Rのパッケージpurrrは、リスト/ベクトルに対して反復処理を行う関数を提供するパッケージである。 これどう呼ぶのが正しいのだろうか。ひとまず「ぴゅるるる」と呼ぼう。

Rでループ処理を行うときは、処理速度が極端に遅いため極力forを使わないようにすべしというのは、ある程度Rを使っている方であればご存知かと思う。 他にはapplyもあるが、個人的にあまり使いやすさを感じたことはない。 purrrは、速度そして可読性の両方を兼ね備えたパッケージである。

基本的な使い方

purrrのメインとなる関数はpurrr::mapである。purrr::map(x, f)として、リスト/ベクトル/行列/データフレームxに対して関数fを繰り返し作用させ、計算結果をリストに格納する。 基本的にはx %>% purrr::map(f)のように、リスト/ベクトルを渡して関数で処理するというパイプラインにするのが一般的かと思う。

例えば、標準正規分布から10個の乱数を取ってきて、それらを5つのリストに格納するコードは以下のようにかける。

library(magrittr)
1:5 %>% purrr::map(rnorm, n = 5)
[[1]]
[1]  1.2644874 -0.1058268  1.1207076  0.9976939  0.1906772

[[2]]
[1] 0.3033158 1.4178929 2.4092054 1.6336716 2.7202050

[[3]]
[1] 2.293005 3.507343 3.918410 2.050489 2.193373

[[4]]
[1] 3.697953 4.015712 5.237348 4.090736 4.329179

[[5]]
[1] 6.763691 4.920744 4.854612 4.336114 5.076095

私の場合、このpurrrを使うシチュエーションは主に2つある。

1. 同じグラフをそれぞれのカテゴリー/グループに対してサクッと作成し、サクッと取り出せるようにしたいとき

1つ目は、同じ種類のグラフ (ヒストグラム、散布図など) を、違うカテゴリー/グループごとに作成したいときである。 ggplot2facet_gridなどを使って、1つの図にまとめる方法もあるが、カテゴリーの数が多すぎるとグラフが小さくなり問題である。 自分の好きなカテゴリーのグラフを手軽に取り出したいときもある。

そこでpurrrの出番である。 irisデータを用いて、アヤメの種類ごとガクの長さ (sepal length) のヒストグラムを作成することを考える。

library(ggplot2)

draw_hist <- function(df, column) {
  g <- ggplot(data = df, aes(x = !!dplyr::enquo(column))) + geom_histogram()
  return(g)
}

hist <- 
  iris %>%
  split(.$Species) %>%
  purrr::map(draw_hist, column = Sepal.Length)

こうすることで、それぞれのアヤメの種類ごとにガクの長さを、例えばhist$setosaのようにして、すぐにsetosaのヒストグラムを呼び出すことができる。

もっというと、irisをlong形式に変換し、各カテゴリー・各変数ごとにsplitしてpurrr::mapしてしまえば、全てのグループと変数の組み合わせについてヒストグラムを作成することができる。

hist_all <- 
  tibble::rowid_to_column(iris, var = "id") %>%
  tidyr::pivot_longer(-c(id, Species), 
                      names_to = "variable", values_to = "value") %>%
  split(list(.$Species, .$variable)) %>%
  purrr::map(draw_hist, column = value)

hist_all$setosa.Petal.Width

データサイズが大きくなったり、分割するカテゴリーが増えたりするとメモリー不足にならないか心配であるが、とりあえずこれまではうまく動いている。

ちなみに上の自作関数の中にある!!dplyr::enquoはnon-standard evaluation (NSE) と呼ばれる黒魔術で、今は無視してほしい。*1

2. 違うデータフレームを続けて結合したいとき

2つ目は、あるデータセットに対して、複数の異なるデータセットを次々と結合させたいときである。 例えば、上で使用したアヤメのデータセットが、アヤメの種類 (Species) 、ガクの長さ (Sepal.Length) 、ガクの幅 (Sepal.Width) 、花弁の長さ (Petal.Length )、花弁の幅 (Petal.Width) の5つに分割されており、これらは全てあるid変数で結合できるとしよう。

df_iris <- tibble::rowid_to_column(iris, var = "id")
df_species <- dplyr::select(df_iris, c(id, Species))
df_sepallength <- dplyr::select(df_iris, c(id, Sepal.Length))
df_sepalwidth <- dplyr::select(df_iris, c(id, Sepal.Width))
df_petallength <- dplyr::select(df_iris, c(id, Petal.Length))
df_petalwidth <- dplyr::select(df_iris, c(id, Petal.Width))

(このコードもpurrrでスマートに書けそうな気がするが、ひとまず力技で押し通した)

これらを元通りに結合するとき、1つのやり方は素直にdplyr::left_joinを繰り返し行うことだろう。

df <- 
  dplyr::left_join(df_species, df_sepallength, by = "id") %>%
  dplyr::left_join(., df_sepalwidth, by = "id") %>%
  dplyr::left_join(., df_petallength, by = "id") %>%
  dplyr::left_join(., df_petalwidth, by = "id")
head(df)
   id Species Sepal.Length Sepal.Width Petal.Length Petal.Width
1  1  setosa          5.1         3.5          1.4         0.2
2  2  setosa          4.9         3.0          1.4         0.2
3  3  setosa          4.7         3.2          1.3         0.2
4  4  setosa          4.6         3.1          1.5         0.2
5  5  setosa          5.0         3.6          1.4         0.2
6  6  setosa          5.4         3.9          1.7         0.4

ただ、コードとしては冗長な感は否めない。どのデータフレームをjoinしたか、混乱することもしばしばある。

私も最近知ったのだが、これはpurrr::reduceを用いるとスッキリと書くことができる。

df <- 
  list(df_species, 
       df_sepallength, 
       df_sepalwidth, 
       df_petallength, 
       df_petalwidth) %>%
  purrr::reduce(dplyr::left_join, by = "id")
head(df)
  id Species Sepal.Length Sepal.Width Petal.Length Petal.Width
1  1  setosa          5.1         3.5          1.4         0.2
2  2  setosa          4.9         3.0          1.4         0.2
3  3  setosa          4.7         3.2          1.3         0.2
4  4  setosa          4.6         3.1          1.5         0.2
5  5  setosa          5.0         3.6          1.4         0.2
6  6  setosa          5.4         3.9          1.7         0.4

繰り返しにはforeachもよい?

forは使うなと冒頭で述べたが、私はちょくちょくforeachを使う。並列処理もforeachで行うことができ、便利である。

標準正規分布から10個の乱数を取ってきて、それらを5つのリストに格納するコードは、foreachを用いて以下のようにも書ける。

library(foreach)
foreach(i = 1:5) %do% {
      rnorm(n = 10)
    }

が、purrrを用いた方が処理速度は速い。それはそうか。

rbenchmark::benchmark(
  use_purrr = {
    1:100 %>% purrr::map(rnorm, n = 1000)
  },
  use_foreach = {
    foreach(i = 1:100) %do% {
      rnorm(n = 1000)
    }
  }
)
         test replications elapsed relative user.self sys.self 
2 use_foreach          100   2.959    3.595     2.935    0.014
1   use_purrr          100   0.823    1.000     0.810    0.011

foreachについてはまた別のポストで書いてみたい。 特に並列処理を行うとき、なぜnamespace::functionと普段から書いておいたほうが良いかがわかる。

*1:いつも雰囲気で何とかしている。

dplyr::acrossを使ってみる

Rでdata frameの特定の列に対してある処理を施すとき、私はdplyr::mutate_atdplyr::mutate_ifをよく使用している。 例えば、本来はnumeric型であるべき列がcharacter型となっているとき、それらの列を全てnumeric型に戻したいときがある。 以下は逆に、dplyr::mutate_ifを用いてnumeric型のものを全てcharacter型にする。

library(magrittr)

iris %>%
  dplyr::mutate_if(.predicate = is.numeric, as.character) %>%
  tibble::as_tibble()

# A tibble: 150 x 5
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
   <chr>        <chr>       <chr>        <chr>       <fct>  
 1 5.1          3.5         1.4          0.2         setosa 
 2 4.9          3           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            3.6         1.4          0.2         setosa 
 6 5.4          3.9         1.7          0.4         setosa 
 7 4.6          3.4         1.4          0.3         setosa 
 8 5            3.4         1.5          0.2         setosa 
 9 4.4          2.9         1.4          0.2         setosa 
10 4.9          3.1         1.5          0.1         setosa 
# … with 140 more rows

これを、dplyr::across関数を使うと、このように書ける。

iris %>%
  dplyr::mutate(dplyr::across(where(is.numeric), as.character)) %>%
  tibble::as_tibble()

# A tibble: 150 x 5
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
   <chr>        <chr>       <chr>        <chr>       <fct>  
 1 5.1          3.5         1.4          0.2         setosa 
 2 4.9          3           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            3.6         1.4          0.2         setosa 
 6 5.4          3.9         1.7          0.4         setosa 
 7 4.6          3.4         1.4          0.3         setosa 
 8 5            3.4         1.5          0.2         setosa 
 9 4.4          2.9         1.4          0.2         setosa 
10 4.9          3.1         1.5          0.1         setosa 
# … with 140 more rows

dplyr::across関数の第1引数として列を指定し、第2引数で処理内容を指定する使い方をするようだ。 今後はdplyr::hogehoge_if, atを使わず、dplyr::acrossを用いるのが推奨されるようだ。

(ちなみに私は、可読性が多少下がるかもしれないがpackage::functionの書き方が好みである。 どのパッケージのどの関数を使っているかを明示的にすること、それによって関数同士のconflictを避けることもできること、並列処理の際に多少楽になることが利点としてあげられる。 このあたりは別途記事にしてみたい)