Sansan Tech Blog

Sansanのものづくりを支えるメンバーの技術やデザイン、プロダクトマネジメントの情報を発信

分析の再現性を担保する工夫

データ分析 R Python

はじめに

技術本部 R&D の小松です。先日、一橋大学の手島健介教授より『経済セミナー』2023年2・3月号をご恵贈いただきました。

www.nippyo.co.jp

手島教授はその中で「米国経済学会データエディター制度の取り組み 再現性向上のためのreplicationチェック」を書かれています。私たちが『経済セミナー』にて「実証研究マネジメントのためのツールキット」の連載時に、手島教授を始めとした研究者の皆さんに草稿を確認いただいたのですが、このトピックはその際に出た議論をまとめられたものです。

そこでは手島教授が体験した、採択された研究論文の再現性チェックのプロセスが事細かに書かれており、興味深いです。現在 AEA P&P のために replication code を準備している私たちにとっても、大変参考になっています。

その中で論文の筆者として行うべきこととして、以下の4点が挙げられています。

  1. AEA が推奨する Social Science Data Editorsのテンプレートを読み、テンプレートに沿って事前に準備しておく。
  2. 秘匿データを使っている場合、どのように説明して、どのように共有するかを考える。
  3. 使っているソフトウェア、パッケージのバージョン、動作環境を確認しておく。
  4. コードが再現性を持つように普段からきちんと準備しておく。

最後の「コードが再現性を持つように普段からきちんと準備しておく」については、トピックから外れるということで軽く触れられるにとどまりました。 しかし、多くの方にとってはどうやって「きちんと準備しておく」かは気になるポイントかと思います。

そこで本稿では、普段の分析で再現性が担保されるように筆者や弊社の他の研究員が工夫している点を簡単に紹介したいと思います。以下、使用する言語は主に R となります。*1

パイプラインツールを用いて分析を行う

再現性を持つように分析を進める方法の1つは、以下の図のように普段からデータの取得 - 加工 - 分析を一貫して回しておくことかと思います。

Source: R for Data Science 2nd Edition.

例えば、データ取得、加工、分析を行うコードをそれぞれ別ファイルにまとめて、それらをたとえば main ファイルですべて回すことが考えられます。

これでも良いと思いますが、不便な点を挙げると

  • 変更を加えていない部分も都度計算され、無駄な待ち時間が発生する
  • 処理時間が長い (数時間 ~ 数日) 処理過程があると、このアプローチは取りづらい
  • そもそも main ファイルに分析コードをまとめられるのは分析プロジェクトが佳境に入った段階であり、探索的にデータを分析しているフェーズではやりづらい

などがあると思います。これら全て、筆者が経験している事象です。

これらを解消する方法として、既存のパイプラインツールを活用することが有り得ます。R であれば、{targets}、 Python であれば Luigi などがあると思います。 Sansanの研究開発部では、Python のプロジェクトでは Luigi のラッパーライブラリである gokart がよく用いられています。

github.com

github.com

github.com

これらライブラリの特徴の1つは、パイプラインの依存関係を記録し、過去すでに行われた計算はスキップし、変更が加えられた部分のみ実行してくれる点です。 R の {targets} については筆者が、Luigi と gokart には弊社研究員の髙橋がそれぞれ過去ブログ記事に書いていますので、そちらも参照下さい。

buildersbox.corp-sansan.com

buildersbox.corp-sansan.com

これらツールを使うと、最初の「変更を加えていない部分も都度計算され、無駄な待ち時間が発生する」というポイントは回避することができます。

次に、「処理時間が長い (数時間 ~ 数日) 処理過程があると、このアプローチは取りづらい」部分についても、パイプラインツールを上手く活用すると効率化することができます。

例として、Johns Hopkins University の Angelo Mele 准教授との研究論文 "Homophily and Community Structure at Scale: An Application to a Large Professional Network" での分析を挙げます。

この研究ではネットワークの各ノードを各クラスター (コミュニティ) に分ける計算過程がありますが、そのクラスターの初期化段階で計算が収束するまで iteration を回す必要があります。

最終的に 20,000 iterations を回すことにしましたが、そこでは 10,000 iteration ごとに分けることにしました。つまり前半に 10,000 iterations を回して、後半 10,001 iteration から再開するイメージです。

...,
  tar_target(
    first_clustering_run,
    lighthergm::hergm(
      g_tokyo_10core ~ edges + triangle + nodematch("h3_7") + nodematch("large_category") + nodematch("occupation_code"),
      n_clusters = 100,
      n_em_step_max = 10000,
      estimate_parameters = FALSE,
      use_infomap_python = TRUE,
      clustering_with_features = FALSE,
      verbose = TRUE,
      infomap_seed = random_seed,
      seeds = c(random_seed)
    )
  ),
  # Apply a few additional iterations of the EM algorithm departing from the checkpoint created by the previous step
  tar_target(
    second_clustering_run,
    lighthergm::hergm(
      first_clustering_run,
      n_em_step_max = 10000,
      clustering_with_features = FALSE
    )
  ),
...

実際に使用しているコードを貼りました。詳細はさておき、iteration を分けて実行しているのがわかるかと思います。

例のように、繰り返し処理を一度に実行するのではなく、分割して実行することで、

  • 段階的に計算結果を確認することができる。例えば、計算が収束しそうか確認したり、計算が意図せず止まった場合でもその止まる前までの結果は参照したりできる。
  • やむを得ず、処理を止めて別の作業をしなければならなくなった場合でも、安心して計算を止めることができる

などのメリットがあると考えられます。寝る前に計算を回し始めて、次の日の朝に計算結果を確認しようとしたら途中で意図しない原因で計算が終わっていたという悲劇を、計算結果を適宜キャッシュしておくことで回避しましょう。

最後に、「そもそも main ファイルに分析コードをまとめられるのは分析プロジェクトが佳境に入った段階であり、探索的にデータを分析しているフェーズではやりづらい 」という点について、{targets} を使っている場合全ての分析をひとつのパイプラインとしてまとめると肥大化しすぎて扱いづらくなるという経験があります。

それを解消するために、{targets} のパイプラインを複数のプロジェクトに分割できる機能の利用を筆者は進めています。 各処理を適切な project の単位で切り出し、別々にパイプラインを実行することで解消できるのではないかと最近トライしているところです。

books.ropensci.org

この辺りはまた知見がまとまり次第、記事としてまとめたいと思います。

分析環境や使用するパッケージを固定する

「使っているソフトウェア、パッケージのバージョン、動作環境を確認しておく」にも関連するポイントですが、分析結果を出すために使用していた環境を、第三者が再現できるようにしておくことは再現性担保のために重要です。

私は R で分析を行うときは、分析環境は Docker を用いて準備し、そこで利用するライブラリは {renv} を用いて管理しています。こうすることで、仮に第三者がある分析を実行したいとなった場合でも、その実行環境を容易に再現することができます。

Docker と {renv} を用いた分析環境の準備については、弊社研究員 Juan Martínez が以下の記事でまとめています。

buildersbox.corp-sansan.com

Docker については、同じくMartínez が以下にまとめていますので、興味のある方はそちらも参照下さい。

justinian336.github.io

データの validation を挟む

実証分析プロジェクトでは、表形式のデータを扱うケースが非常に多いと思います。そのデータが意図した通りのデータになっているかを都度確認するプロセスを挟んでおくと、思わぬバグを防ぐことができます。

Python で表形式のデータを扱うときは pandas の DataFrame object を扱うことが多いと思います。そのデータを validation する上で使うpandera があります。

github.com

弊社の研究開発部でもデータフレームを扱うプロジェクトではよく使用しています。 例えば、あるアンケート回答者のデータフレームを扱うとき、そのデータフレームの「型」 (スキーマ (schema) とも言ったりします) を以下のように表現することができます。

import pandera as pa
from pandera.typing import Series


class RespondentSchema(pa.SchemaModel):
    id: Series[str] = pa.Field(nullable=False, unique=True, description="Respondent ID for each record")
    prefecture: Series[str] = pa.Field(nullable=False,  description="Prefecthre where each respondent lives")
    age: Series[int] = pa.Field(nullable=True, ge=0, description="Age of each respondent")

こうすることで、そのデータフレームの各列の型 (string なのか ingeger なのか、ユニークなのかどうか、欠損を許すかどうか、値の範囲はどこまでか…) を表現することができます。 表現するだけでなく、validate メソッドを使うことで

RespondentSchema.validate(df)

のように、扱うデータフレームが定義した型に違反していないかを確認することができます。

pandera については、弊社研究員の前嶋が過去に記事にしているので、よろしければそちらもご覧ください。

buildersbox.corp-sansan.com

R にもこういったデータを validation するためのライブラリはいくつか存在します。最近使用を検討しているのは {pointblank} というライブラリです。

github.com

上の pandera で書いたようなスキーマは

validate_respondent_data <- function(data) {
  data |>
    pointblank::col_exists(c("id", "prefecture", "age")) |>
    pointblank::col_is_character(c("id", "prefecture")) |>
    pointblank::col_is_integer("age") |>
    pointblank::rows_complete(c("id", "prefecture")) |>
    pointblank::rows_distinct(c("id")) |>
    pointblank::col_vals_gte("age", 0, na_pass = TRUE)
}

と書けたりします。ただし pandera と比較すると少し冗長ですね。*2

このデータフレームの validation が特に効くケースとして、データ同士を id を用いて結合するものがあるでしょう。　

df1 <- tibble::tibble(
  id = c("A", "B"),
  feature_x = c(1, 2)
)

df2 <- tibble::tibble(
  id = c("A", "A", "B"),
  feature_y = c(334, 334, 256)
)

dplyr::left_join(df1, df2, by = "id")

ここでは df1 に df2 を結合しようとしています。分析者が期待する出力は

# A tibble: 3 × 2
  id    feature_x feature_y
  <chr>     <dbl>     <dbl>
1 A             1       334
3 B             2       256

ですが、df2 に分析者が意図していない重複が存在しているため、実際には

# A tibble: 3 × 3
  id    feature_x feature_y
  <chr>     <dbl>     <dbl>
1 A             1       334
2 A             1       334
3 B             2       256

という A についてレコードが重複してしまったデータフレームが出力されています。これはよくあるミスではないでしょうか。これは、

validate_df2 <- function(data) {
  data |>
    pointblank::rows_distinct(c("id"))
}

のような validator を定義し、df2 を事前に validation して id に重複が無いかを確認しておけば防ぐことができます。 加えて、結合したあとのデータフレームについても、id が重複していないかなど意図した通りにデータフレームが生成されているか validation を挟んでおくと良いでしょう。

これらの工夫は、再現性の担保というよりかは、分析結果の「正しさ」を担保するために役に立つ機構であると言えそうです。

処理を関数化し、テストする

分析を進める過程で、似たような処理を何度も実行することが増えてくるかと思います。例えば、

df <- tibble::tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)

df$a <- (df$a - min(df$a, na.rm = TRUE)) / 
  (max(df$a, na.rm = TRUE) - min(df$a, na.rm = TRUE))
df$b <- (df$b - min(df$b, na.rm = TRUE)) / 
  (max(df$b, na.rm = TRUE) - min(df$a, na.rm = TRUE))
df$c <- (df$c - min(df$c, na.rm = TRUE)) / 
  (max(df$c, na.rm = TRUE) - min(df$c, na.rm = TRUE))
df$d <- (df$d - min(df$d, na.rm = TRUE)) / 
  (max(df$d, na.rm = TRUE) - min(df$d, na.rm = TRUE))

という処理は、

normalize_between_zero_one <- function(x) {
  ( x - max(x, na.rm = TRUE) ) / ( max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
}

という関数を定義して

df <- df |>
  dplyr::mutate(
    dplyr::across(c(a, b, c, d), normalize_between_zero_one)
  )

としてあげた方が読みやすいし、無駄な copy & paste が発生しないのでコピペ間違いによるバグも防ぐことができそうです。 自分はそんなことをしていないと思っていても、自分の分析コードを読み返すとこういう箇所はしばしばあったりします。筆者もかつては最初のような重複が沢山あるコードを書いていましたし、今も気を抜くとすぐに書いてしまいます。

関数化するメリットは重複する処理の効率化の他に、その挙動を確認するために単体テスト (unit testing) を書くことが簡単になることです。 単体テストとは、例えば

make_from_start_to_end_yearmonth <- function(start_yearmonth, end_yearmonth) {
  seq(
    as.Date(paste0(start_yearmonth, "-01")),
    as.Date(paste0(end_yearmonth, "-01")),
    by = "months"
  ) %>%
    format("%Y-%m") %>%
    as.character()
}

という関数があったときに、この関数が意図した通りに動くかを検証するテストのことを言います。R では {testthat} というパッケージを使って

test_that("test_make_from_start_to_end_yearmonth", {
  output_data <- make_from_start_to_end_yearmonth("2022-10", "2023-02")
  expected_data <- c("2022-10", "2022-11", "2022-12", "2023-01", "2022-02")
  testthat::expect_equal(output_data, expected_data)
})

のように書くことができます。

この単体テストという考え方は、研究者の方にはあまり馴染みがないかもしれません。私も大学院生のときは知りませんでしたし、今の会社に入ってその存在を知ったときは、 なぜそんなテストを書く必要があるのかとその存在意義に懐疑的でした。しかしテストを書く習慣をつけると、良いことがいくつかあります。

1つ目は、バグを減らすことができる点です。これはテストを書く目的が関数が意図した通り動くかを確認するためですから、その通りですね。テストを通じてバグに対して自覚的になることは、その早期発見につながります。

2つ目は、単体テストは定義した関数がどのような挙動をするのかというドキュメンテーションも兼ねる点です。 関数の定義を見ても、入力が何で出力が何かがわかりにくいこともありますが、単体テストが書かれていると、その input-output の関係がクリアになります。それは上の make_from_start_to_end_yearmonth() の例を見ていただくとよく分かるかと思います。

最後は、より良いコードを書く手助けになることです。ある関数の単体テストを書こうとしたとき、テストが書きづらいと思うことがあります。それは、もとの関数があまり上手く書けていない可能性が高いです。

例を挙げます。例えば以下のような適当な関数を定義します。*3

clean_data <- function(data, col_str, col_group, col_val) {
  data |>
    dplyr::mutate({{ col_str }} := stringr::str_replace_all({{ col_str }}, "[[:blank:]]", "")) |>
    dplyr::group_by({{ col_group }}) |>
    dplyr::arrange({{ col_value }}) |>
    dplyr::slice(1) |>
    dplyr::ungroup()
}

この関数は

  1. col_str に whitespace が含まれている場合は、それを全て取り除く
  2. 各col_group ごとに、最大の col_value を持つ行を残す

以上2つの処理を行っています。複数の作業を1つの関数で一度に行っているため、whitespace の処理に関するテストと、行を絞る処理のテストを別々に書くことが難しくなっています。 これがテストが「書きづらい」例かと思います。一般化は難しいですが、1つの関数に複数の処理が行われると書きづらいケースが多いように思います。

このときは、上で書いた2つの処理を別々にすることを考えます。つまり、

remove_blank_from_string_column <- function(data, col_str) {
  data |>
    dplyr::mutate({{ col_str }} := stringr::str_replace_all({{ col_str }}, "[[:blank:]]", ""))
}

keep_row_with_largest_value_for_each_group <- function(data, col_group, col_value) {
  data |>
    dplyr::group_by({{ col_group }}) |>
    dplyr::arrange(desc({{ col_value }})) |>
    dplyr::slice(1) |>
    dplyr::ungroup()
}

と分けてあげることで、

test_that("test_remove_blank_from_string_column", {
  test_data <- tibble::tibble(x = c("A", "B ", "C  ", "D　", "E E "))
  output_data <- remove_blank_from_string_column(test_data, x)
  expected_data <- tibble::tibble(
    x = c("A", "B", "C", "D", "EE")
  )
  testthat::expect_equal(output_data, expected_data)
})

test_that("test_remove_blank_from_string_column", {
  test_data <- tibble::tribble(
    ~y,  ~z,
    "a", 1,
    "b", 1,
    "c", 1,
    "a", 10,
    "b", 1
  )
  output_data <- keep_row_with_largest_value_for_each_group(test_data, y, z)
  expected_data <- test_data <- tibble::tribble(
    ~y,  ~z,
    "a", 10,
    "b", 1,
    "c", 1
  )
  testthat::expect_equal(output_data, expected_data)
})

のようにテストが書きやすくなりました。このように、1つの関数にあまり多くの処理を持たせすぎない、もっと言えば1つの関数が持つべき役割が1つだけになるようにしてあげると、見通しがよくなります。 これは 単一責任の原則 (Single responsibility principle) とも呼ばれます。テストを書くことを前提にすると、良いコードを書くことにつながりやすいです。

この上で、もとの clean_data() と同じ挙動をする関数を作りたい場合は、いろいろ方法はあると思いますが、例えば関数型プログラミングの考え方を使えば

clean_data <- purrr::compose(
  !!!list(
    purrr::partial(remove_blank_from_string_column, col_str = x),
    purrr::partial(keep_row_with_largest_value_for_each_group, col_group = y, col_value = z)
  ), 
  .dir = "forward"
)

と書くことができるでしょう。*4

単体テストの実行は、R ではパッケージ化することでよりやりやすくなります。詳しくは R Packages (2e) をご覧になると良いでしょう。*5

以上述べた関数化と単体テストの考え方も、どちらかと言えば分析の「正しさ」を担保する工夫と言えるかもしれません。

最後に

以上簡単ですが、分析の再現性を担保するための工夫を紹介しました。まだまだ模索している段階であり、他にも工夫の余地はたくさんあるに違いありません。そういったアイデアがあればシェアいただけると嬉しく思います。

この記事を書く過程で、弊社研究員の前嶋の以下の Tweet を思い出しました。

私が紹介した工夫は、ソフトウェア・エンジニアリング界隈では馴染みのある話がほとんどかと思いますが、(経済学の) 実証研究に従事している方々にとっては慣れないものもあるのかなと思います。

では研究者がソフトウェア・エンジニアリングを学ぶべきかと言えば、学ぶに越したことはありませんけれども、それよりも前嶋の言うように役割分担して各々得意な領域で貢献できるような世界であれば、皆幸せになれるのかもしれません。

筆者は研究者としてもソフトウェア・エンジニアとしても大したことはありませんが、それら両方にそこそこ従事している身として上手に橋渡しできる役回りとして仕事ができないかなと思ったりしています。もしご興味あればお声がけ下さいませ。

*1:Stata ユーザーの皆さん、ごめんなさい。

*2:類似する R のパッケージは他には validate があります。ただ、pandera に比べるとスキーマが書きにくい印象があり、色々試行錯誤しています。より良いパッケージをご存知でしたら教えて下さい。

*3:ここで使用している {{}} は curly-curly operator と呼ばれるもので、tidyverse 系の関数を自作関数で用いるときに使うと便利です。 https://www.tidyverse.org/blog/2019/06/rlang-0-4-0/#a-simpler-interpolation-pattern-with-

*4:例として不適切だったかもしれません。

*5:研究では、使う関数をパッケージ化しながらパイプラインを組んでいます。

© Sansan, Inc.