はじめに

twitterで以下のツイートを見かけました。

R :: mixgb 📦https://t.co/TuD7bEU9hT

mixgb は XGBoost を使用した多重代入法の R 📦

XGBoost による予測値に対して，predictive mean matching (PMM) で補完をしており，https://t.co/BYkjF6eg0E
にあるような複数の方法が実装されている（Algo. 1）

🔖 arXiv: https://t.co/wmwVcMyeB7 pic.twitter.com/gelg2hJYF5
— Maxwell! 🌃 🌻 (@Maxwell_110) 2022年8月4日

リポジトリはこちらです。

github.com

機械学習による多重代入というものが気になったので、論文を読んで簡単にまとめます。

概要

Rubin (1978)によって多重代入法(Multiple imputation:MI)が提案された。そしてMIを実施するために、多くのフレームワークが開発されてきている。

従来のMI実装の主な欠点の1つは、変数間の複雑な関係を捉えることができないこと。

これを解決するためにツリーベースのアルゴリズムが提案されている。 Stekhoven and Bühlmann (2012)は、ランダムフォレストに基づく欠損値代入のノンパラメトリック手法を提案し、missForestというパッケージで実装している。

もう一つの欠点は、大規模なデータセットに対して計算コストが大きいこと。

こちらに関しても近年の技術の進歩により、大規模なデータセットの収集と分析が可能になったが、ツリーベースのモデルを利用する手法は小規模なデータセットに適しており、中規模のデータセットへの適用が大変。

またMLベースの代入の欠点としては、単一値の代入しかできず欠測値に関する不確実性を加味できない。

そこで今回の論文では、大規模なデータセットに対しても高速な処理を行うことができ、欠損値の不確実性を加味したXGBoostフレームワークによる多重代入（mixgb）を提案している。

xgboostによる多重代入法

通常の機械学習による代入は欠損データの不確実性を十分に考慮することなく、欠損値の予測精度に焦点を当てている。

その結果、代入の変動性を過小評価することになる(らしい)。

この問題に対処するため、ブートストラップ法と予測平均マッチング（PMM）を用いたXGBoostによる代入を検討している。

大まかな流れは以下の通り

各特徴量の欠損割合が昇順になるように、特徴量をソートする。

全体サンプル $Y$ からブートストラップによるサブサンプル $Y^*$ を取得する
代入対象の特徴量以外の特徴量 $Y _ { -i }^ {*obs }$ で、代入対象の特徴量 $Y _ { i } ^ { *obs }$ を予測するXGBoostモデル $f^*$ を作成する。このモデルを用いて欠測している特徴量 $\widetilde{Y_i^{mis}}$ を予測する。
PMMの実施

　type1の場合

観測値全体を用い代入対象の特徴量 $i$ を予測するXGBoostモデル $f$ を作成
このモデル $f$ を用いて観測値の推定量 $\widehat{ Y_i^{obs}}$ を算出
$\widehat{ Y_i^{obs}}$ と $\widetilde{Y_i^{mis}}$ を用いてPMMを行う

　type2の場合

モデル $f^*$ を用いて観測値の推定量 $\widetilde{ Y_i^{obs}}$ を算出
$\widehat{ Y_i^{obs}}$ と $\widetilde{Y_i^{mis}}$ を用いてPMMを行う

　4.1-2の処理を各特徴量に対して順に適用して欠測値の代入していく。

そして、この1~4の代入操作をM回繰り返す

シミュレーション

シミュレーションを行い複数の手法におけるバイアスの大きさなどを確認している。

利用したデータは以下の式に従い生成している。

それぞれR=1000サンプルを生成し、これらの変数の一部に20%の割合でMAR, MCARで欠測させている。

確認しているのはBiasとVarianceである。

論文で掲載されている結果では、BiasとVariance共にMICEやmisForestよりも小さい傾向にあると考えられるが、変数と設定次第で変化し、絶対的に良い代入ができているというわけではなさそうだ。

Rによる実装

mixbgのサンプルを実行して挙動を確認します。

はじめにインストールと読み込み。

install.packages("mixgb")
library(mixgb)

mixgbパッケージにはnhanes3_newbornと呼ばれる1988-1994年の新生児のデータセットが含まれている。

性別や体重、頭囲、肩甲骨部や上腕三頭筋の皮下組織の厚さなどの情報が入っている。

str(nhanes3_newborn)

##  tibble [2,107 × 16] (S3: tbl_df/tbl/data.frame)
##   $ HSHSIZER: int [1:2107] 4 3 5 4 4 3 5 3 3 3 ...
##   $ HSAGEIR : int [1:2107] 2 5 10 10 8 3 10 7 2 7 ...
##   $ HSSEX   : Factor w/ 2 levels "1","2": 2 1 2 2 1 1 2 2 2 1 ...
##   $ DMARACER: Factor w/ 3 levels "1","2","3": 1 1 2 1 1 1 2 1 2 2 ...
##   $ DMAETHNR: Factor w/ 3 levels "1","2","3": 3 1 3 3 3 3 3 3 3 3 ...
##   $ DMARETHN: Factor w/ 4 levels "1","2","3","4": 1 3 2 1 1 1 2 1 2 2 ...
##   $ BMPHEAD : num [1:2107] 39.3 45.4 43.9 45.8 44.9 42.2 45.8 NA 40.2 44.5 ...
##    ..- attr(*, "label")= chr "Head circumference (cm)"
##   $ BMPRECUM: num [1:2107] 59.5 69.2 69.8 73.8 69 61.7 74.8 NA 64.5 70.2 ...
##    ..- attr(*, "label")= chr "Recumbent length (cm)"
##   $ BMPSB1  : num [1:2107] 8.2 13 6 8 8.2 9.4 5.2 NA 7 5.9 ...
##    ..- attr(*, "label")= chr "First subscapular skinfold (mm)"
##   $ BMPSB2  : num [1:2107] 8 13 5.6 10 7.8 8.4 5.2 NA 7 5.4 ...
##    ..- attr(*, "label")= chr "Second subscapular skinfold (mm)"
##   $ BMPTR1  : num [1:2107] 9 15.6 7 16.4 9.8 9.6 5.8 NA 11 6.8 ...
##    ..- attr(*, "label")= chr "First triceps skinfold (mm)"
##   $ BMPTR2  : num [1:2107] 9.4 14 8.2 12 8.8 8.2 6.6 NA 10.9 7.6 ...
##    ..- attr(*, "label")= chr "Second triceps skinfold (mm)"
##   $ BMPWT   : num [1:2107] 6.35 9.45 7.15 10.7 9.35 7.15 8.35 NA 7.35 8.65 ...
##    ..- attr(*, "label")= chr "Weight (kg)"
##   $ DMPPIR  : num [1:2107] 3.186 1.269 0.416 2.063 1.464 ...
##    ..- attr(*, "label")= chr "Poverty income ratio"
##   $ HFF1    : Factor w/ 2 levels "1","2": 2 2 1 1 1 2 2 1 2 1 ...
##   $ HYD1    : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 3 1 1 1 1 1 1 2 1 ...

VIMパッケージを利用して欠測状況を確認する。　

# 欠損状況の確認
library(VIM)
vim.aggr <- aggr(nhanes3_newborn,
                 col = c('white','red'),
                 numbers = TRUE,
                 sortVars = FALSE,
                 prop = TRUE,
                 labels = names(nhanes3_newborn),
                 cex.axis = .8,
                 gap = 5)

一部欠測状況に偏りがあるように見受けられる。

代入処理を行う前の以下のようにデータクリーニングを行う必要がある。

データはデータフレームで与えられる
IDを削除する必要がある
欠損値はNaN
Infや-Infは含めない
Char型変数は、Factor型に変換する
Factor型の変数は、少なくとも 2 つのレベルが必要

このデータに対してmixgb関数を用いてXGBoostによる代入を行う。

# 代入の実施
imputed_data <- mixgb(data = nhanes3_newborn, m = 5)

mixgb関数の引数には次のように設定できる。

m：出力データセットの数
maxit：代入の反復回数(デフォルト:1)
ordinalAsInteger：順序変数を整数に変換するかどうか
bootstrap：ブートストラップを使用するかどうか
pmm.type：予測平均マッチング設定

そのほかにもXGBoostモデルのハイパーパラメータの設定を設定することができる。

代入結果の状態を確認する。いくつかの描画関数が用意されている。

はじめにヒストグラムの確認

plot_hist(imputation.list = imputed_data,
          var.name = "DMPPIR",
          original.data = nhanes3_newborn)

次に２変数(連続値×連続値)の分布の確認

plot_2num(imputation.list = imputed_data, 
          var.x = "BMPSB2",
          var.y = "BMPSB1", 
          original.data = nhanes3_newborn)

そして２変数(連続値×離散値)の分布の確認

plot_1num1fac(imputation.list = imputed_data, 
              var.num = "BMPSB2",
              var.fac = "HSSEX", 
              original.data = nhanes3_newborn)