mlrパッケージを使った欠損値代入

以下の文章はmlrパッケージのチュートリアルの一部である[Imputation - mlr tutorial](https://mlr-org.github.io/mlr-tutorial/devel/html/impute/index.html)を訳したものです。可読性のためにコードと文章は若干変更しています。

欠損値に対する代入
==================

欠損値に対してmlrパッケージは、[imputations function | R Documentation](https://www.rdocumentation.org/packages/mlr/versions/2.12.1/topics/imputations)にリストアップしている複数の代入法をサポートしている。この中には、定数(固定値、平均値、中央値、最頻値など)で代入する標準的な方法と、乱数(対象の特徴量の経験分布または特定の分布族からの抽出)を用いる方法の両方が含まれている。さらに、欠損値を他の特徴量に基づく教師あり学習からの予測値で置き換える方法もmlrには含まれている。

もし望みの代入方法が存在しなければ、これを拡張することも簡単にできる(詳しくは[Create Imputation Methods - mlr tutorial](https://mlr-org.github.io/mlr-tutorial/devel/html/create_imputation/index.html)を参照のこと)。

mlrに含まれている学習アルゴリズムの中には、欠損値を適当な方法(つまり、単に削除する以外の方法で)扱えるものも含まれているという点にも注意してほしい。そのような学習器は"missings"プロパティを持っている。これは、`listLearners`関数を用いると確認することができる。

``` r
listLearners("regr", properties = "missings")[c("class", "package")]
```

    ##              class      package
    ## 1 regr.bartMachine  bartMachine
    ## 2  regr.blackboost mboost,party
    ## 3     regr.cforest        party
    ## 4       regr.ctree        party
    ## 5      regr.cubist       Cubist
    ## 6 regr.featureless          mlr
    ## ... (#rows: 14, #cols: 2)

詳細については[Integrated Learners - mlr tutorial](https://mlr-org.github.io/mlr-tutorial/devel/html/integrated_learners/index.html)を確認してもらいたい。

`impute`関数と`reimpute`関数
------------

代入は`impute`関数により実行される。代入方法は、特徴量ごとに、あるいは特徴量のクラス(数値型なのか因子型なのか)ごとに指定することができる。さらに、特徴量あるいはクラスのいずれを指定する場合でも、欠損値を示すダミー変数の生成を指示できる。これにより、欠損値のデータの傾向を特定することや、後の分析で代入された値と本来の観測値を分けて取り扱うことが可能となる。

まず、例として`airquality`データセットを見てみよう(cf. [airquality function | R Documentation](https://www.rdocumentation.org/packages/datasets/versions/3.5.0/topics/airquality))。

``` r
summary(airquality)
```

    ##      Ozone           Solar.R           Wind             Temp      
    ##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
    ##  1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
    ##  Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
    ##  Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
    ##  3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
    ##  Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
    ##  NA's   :37       NA's   :7                                       
    ##      Month            Day      
    ##  Min.   :5.000   Min.   : 1.0  
    ##  1st Qu.:6.000   1st Qu.: 8.0  
    ##  Median :7.000   Median :16.0  
    ##  Mean   :6.993   Mean   :15.8  
    ##  3rd Qu.:8.000   3rd Qu.:23.0  
    ##  Max.   :9.000   Max.   :31.0  
    ## 

`Ozone`には37、`Solar`には7つの欠損値があることが分かる。今回、例のために`Wind`に人工的な欠損値を追加し、さらに`Wind`を因子型変数に変換する。

``` r
set.seed(123)
airq <- airquality
ind <- sample(nrow(airq), 10)
airq$Wind[ind] <-  NA
airq$Wind <- cut(airq$Wind, c(0, 8, 16, 24))
summary(airq)
```

    ##      Ozone           Solar.R           Wind         Temp      
    ##  Min.   :  1.00   Min.   :  7.0   (0,8]  :50   Min.   :56.00  
    ##  1st Qu.: 18.00   1st Qu.:115.8   (8,16] :86   1st Qu.:72.00  
    ##  Median : 31.50   Median :205.0   (16,24]: 7   Median :79.00  
    ##  Mean   : 42.13   Mean   :185.9   NA's   :10   Mean   :77.88  
    ##  3rd Qu.: 63.25   3rd Qu.:258.8                3rd Qu.:85.00  
    ##  Max.   :168.00   Max.   :334.0                Max.   :97.00  
    ##  NA's   :37       NA's   :7                                   
    ##      Month            Day      
    ##  Min.   :5.000   Min.   : 1.0  
    ##  1st Qu.:6.000   1st Qu.: 8.0  
    ##  Median :7.000   Median :16.0  
    ##  Mean   :6.993   Mean   :15.8  
    ##  3rd Qu.:8.000   3rd Qu.:23.0  
    ##  Max.   :9.000   Max.   :31.0  
    ## 

まず、全ての整数型特徴量の欠損値を平均値で、因子型特徴量の欠損値を最頻値で代入し、全ての整数型特徴量に対してダミー変数を生成する方法を示す。

``` r
imp = impute(
  airq, # データセットの指定
  classes = list( # クラスごとに代入方法を指定
    integer = imputeMean(), 
    factor = imputeMode()
  ), 
  dummy.classes = "integer" # ダミー変数を作成するクラスの指定
)
```

`impute`の返り値はリストであり、`$data`スロットに代入済みのデータセットが含まれる。デフォルトでは、ダミー変数は`TRUE`と`FALSE`の2水準を持つ因子型である。これはオプションで0と1の変数にすることもできる。

``` r
head(imp$data)
```

    ##      Ozone  Solar.R   Wind Temp Month Day Ozone.dummy Solar.R.dummy
    ## 1 41.00000 190.0000  (0,8]   67     5   1       FALSE         FALSE
    ## 2 36.00000 118.0000  (0,8]   72     5   2       FALSE         FALSE
    ## 3 12.00000 149.0000 (8,16]   74     5   3       FALSE         FALSE
    ## 4 18.00000 313.0000 (8,16]   62     5   4       FALSE         FALSE
    ## 5 42.12931 185.9315 (8,16]   56     5   5        TRUE          TRUE
    ## 6 28.00000 185.9315 (8,16]   66     5   6       FALSE          TRUE

`$desc`スロットにはImputationDescオブジェクトが格納されている。このオブジェクトの中には代入に関する全ての情報が保存されている。今回の例では、代入のために計算された平均値や最頻値が含まれている(訳注: `imp$desc$impute`の中に計算結果が入っているが表示が長いので略)。

``` r
imp$desc
```

    ## Imputation description
    ## Target: 
    ## Features: 6; Imputed: 6
    ## impute.new.levels: TRUE
    ## recode.factor.levels: TRUE
    ## dummy.type: factor

上記Imputation descriptionの中には、目的変数の名前(今の例では存在しない)、特徴量の数と代入された特徴量の数が表示される。なお、ここでいう**代入された特徴量の数**というのは、実際にNAを含んでいて代入が行われた特徴量の数**という意味ではない**ことに注意してほしい。この数は、代入方法が指定された特徴量の数である。今回は、5つの整数型特徴量と1つの因子型特徴量に対して代入方法を指定しているので、6という値になっている。`dummy.type`は、作成されたダミー変数が因子型であることを示している。`impute.factor.levels`と`recode.factor.levels`についての詳細はimpute関数のヘルプ([impute function | R Documentation](https://www.rdocumentation.org/packages/mlr/versions/2.12.1/topics/impute))を参照してもらいたい。

次に、目的変数を含む場合の例を見てみよう。airqualityデータを使った学習タスクの例として、オゾン濃度を他の気象要素から予測する、というものを考えてみる。この目的のためには`Day`と`Month`列は不要なので、まずはこれを除外する。

``` r
airq <- subset(airq, select = 1:4)
```

最初の100の観測値を訓練セット、残りをテストセットとする。

``` r
airq.train <- airq[1:100, ]
airq.test <- airq[-c(1:100), ]
```

教師あり学習の問題においては、impute関数を呼び出す際に目的変数を指定する必要がある。これにより、代入とダミー変数の生成が目的変数に適用されなくなり、(他の特徴量に対する)代入のために目的変数が使われることもなくなる。

先程は代入する特徴量をクラス単位で指定したが、今回は特徴量毎に代入方法を指定してみよう。

`Solar.R`の欠損値は、欠損していない観測値に基づく経験分布からの乱数で代入する。

`inputeLearner`関数は、代入のためにmlrに統合された全ての教師あり学習アルゴリズムを使用することができる。学習器のタイプ(`regr`なのか`classif`なのか)は、特徴量のクラスに応じて指定しなければならない。今回、`Wind`の欠損値は分類木(`rpart`)による予測結果で代入する。このとき、デフォルトでは目的変数と`Wind`以外の利用可能な全ての特徴量が分類木構築のために使われる(今回の例では`Solar.R`と`Temp`だ)。使用する特徴量は任意に指定することもできる。なお、`rpart`は欠損値があっても動作するので、`Solar.R`のNAは問題にはならない。

``` r
imp <- impute(
  airq.train, 
  target = "Ozone", 
  cols = list(
    Solar.R = imputeHist(),
    Wind = imputeLearner("classif.rpart")
  ),
  dummy.cols = c("Solar.R", "Wind")
)
summary(imp$data)
```

    ##      Ozone           Solar.R           Wind         Temp      
    ##  Min.   :  1.00   Min.   :  7.0   (0,8]  :33   Min.   :56.00  
    ##  1st Qu.: 16.00   1st Qu.:100.5   (8,16] :61   1st Qu.:69.00  
    ##  Median : 34.00   Median :223.0   (16,24]: 6   Median :79.50  
    ##  Mean   : 41.59   Mean   :194.0                Mean   :76.87  
    ##  3rd Qu.: 63.00   3rd Qu.:274.2                3rd Qu.:84.00  
    ##  Max.   :135.00   Max.   :334.0                Max.   :93.00  
    ##  NA's   :31                                                   
    ##  Solar.R.dummy Wind.dummy
    ##  FALSE:93      FALSE:94  
    ##  TRUE : 7      TRUE : 6  
    ##                          
    ##                          
    ##                          
    ##                          
    ## 

``` r
imp$desc
```

    ## Imputation description
    ## Target: Ozone
    ## Features: 3; Imputed: 2
    ## impute.new.levels: TRUE
    ## recode.factor.levels: TRUE
    ## dummy.type: factor

ImputationDescオブジェクトは、`reimpute`関数に引数として与えることで、テストセットに訓練セットと同様の方法で代入を行うことができる。

``` r
airq.test.imp <- reimpute(airq.test, imp$desc)
head(airq.test.imp)
```

    ##   Ozone Solar.R   Wind Temp Solar.R.dummy Wind.dummy
    ## 1   110     207  (0,8]   90         FALSE      FALSE
    ## 2    NA     222 (8,16]   92         FALSE      FALSE
    ## 3    NA     137 (8,16]   86         FALSE      FALSE
    ## 4    44     192 (8,16]   86         FALSE      FALSE
    ## 5    28     273 (8,16]   82         FALSE      FALSE
    ## 6    65     157 (8,16]   80         FALSE      FALSE

いくつかのリサンプリング手法によって機械学習の手法を評価したいという場合には、impute/reimpute関数が、訓練と予測の前に自動的に呼び出されるのが望ましいと考えるだろう。これは、代入ラッパーを作成することで実現できる。

学習器と代入の融合
------------------

`makeImputeWrapper`関数を使うことで、学習器と代入を組み合わせることができる。この関数の引数は基本的には`impute`関数と同様である。先ほどと同様に、`Solar.R`は経験分布に基づく乱数で、`Wind`は分類木による予測値で代入し、さらに両方の変数について代入した値であるかどうかを示すダミー変数を生成しよう。

``` r
lrn = makeImputeWrapper(
  "regr.lm",
  cols = list(
    Solar.R = imputeHist(),
    Wind = imputeLearner("classif.rpart")
  ),
  dummy.cols = c("Solar.R", "Wind")
)
lrn
```

    ## Learner regr.lm.imputed from package stats
    ## Type: regr
    ## Name: ; Short name: 
    ## Class: ImputeWrapper
    ## Properties: numerics,factors,se,weights,missings
    ## Predict-Type: response
    ## Hyperparameters:

学習器を訓練する前には、`impute`関数が呼び出されて訓練セットに対する代入が行われる。また、予測を行う前には`reimpute`関数が呼び出され、訓練段階で作成されたImputationDescオブジェクトを使ってテストセットに対する代入が行われる。

さて、再びオゾン濃度を気象要素から予測することを考えてみよう。タスクを作成するためには、目的変数に欠損値があるレコードを削除しなければならない。

``` r
airq <- subset(airq, subset = !is.na(airq$Ozone))
task <- makeRegrTask(data = airq, target = "Ozone")
```

続いて、3分割クロスバリデーションによって平均二乗誤差を計算する。

``` r
rdesc <- makeResampleDesc("CV", iter = 3)
r <- resample(lrn, task, resampling = rdesc, show.info = FALSE, models = TRUE)
r$aggr
```

    ## mse.test.mean 
    ##      591.4208

``` r
lapply(r$models, getLearnerModel, more.unwrap = TRUE)
```

    ## [[1]]
    ## 
    ## Call:
    ## stats::lm(formula = f, data = d)
    ## 
    ## Coefficients:
    ##       (Intercept)            Solar.R         Wind(8,16]  
    ##        -108.30884            0.08772          -20.67284  
    ##       Wind(16,24]               Temp  Solar.R.dummyTRUE  
    ##          -5.34114            1.88699            3.92753  
    ##    Wind.dummyTRUE  
    ##         -13.31178  
    ## 
    ## 
    ## [[2]]
    ## 
    ## Call:
    ## stats::lm(formula = f, data = d)
    ## 
    ## Coefficients:
    ##       (Intercept)            Solar.R         Wind(8,16]  
    ##         -79.88831            0.02543          -28.12074  
    ##       Wind(16,24]               Temp  Solar.R.dummyTRUE  
    ##         -25.43269            1.74266          -22.06698  
    ##    Wind.dummyTRUE  
    ##          10.07808  
    ## 
    ## 
    ## [[3]]
    ## 
    ## Call:
    ## stats::lm(formula = f, data = d)
    ## 
    ## Coefficients:
    ##       (Intercept)            Solar.R         Wind(8,16]  
    ##         -95.86501            0.06464          -12.92908  
    ##       Wind(16,24]               Temp  Solar.R.dummyTRUE  
    ##         -12.58382            1.71833          -12.36961  
    ##    Wind.dummyTRUE  
    ##           4.22593

この他に、`makePreprocWrapperCaret`関数をつかって学習器と代入を融合する手段もある。この関数は`caret`パッケージの`preProcess`関数へのインターフェースとなっている。`preProcess`関数がサポートするのは数値型の特徴量のみであり、k最近傍点法、bagged tree、中央値による代入が可能である。