1. エンジン名定義

今回は線型回帰（linear_reg）にロバスト線型回帰（rlm）を追加したいので、 以下のように引数を指定します。

parsnip::set_model_engine(
  model = "linear_reg",
  mode = "regression",
  eng = "rlm"
)

mode引数は回帰の目的が連続値（"regression"）かクラス（"classification"）かを指定します。

2. 依存関係定義

rlm関数はMASSパッケージが提供しているので、以下のように定義します。

parsnip::set_dependency(
  model = "linear_reg",
  eng = "rlm",
  pkg = "MASS"
)

3. 学習方法定義

これまでと同様にmodel、eng、modeの引数を指定します。

特殊なのはvalue引数でしょう。 ここには、interface, protect, func, defaultsの4つの名前を持つリストを与えます¹。

• interface
  • オリジナル実装の学習用関数が、データをどのように受け取るか指定します。
    • formula: formula引数にモデル式を、data引数にデータフレームを入力する場合
    • data.frame: x引数に特徴量のデータフレームを、y引数にベクトルを指定する場合
    • matrix: x引数に計画行列を、yにベクトルを指定する場合
  • いずれのinterfaceも引数の名前と順序が固定されている点に注意してください。 たとえば引数がyとXの順であれば、いずれのインターフェースにも摘要しません。 その場合、オリジナルの関数をいずれかのインターフェースに適合する関数でラップし、パッケージ化が必要です。
• protect
  • ユーザーが指定不能な引数
    • interfaceが利用する引数
    • weightsなど重みを指定する引数
      • 現段階では重み付けに対応していないため
      • https://github.com/tidymodels/parsnip/issues/136
    • その他
  • 例: c("formula", "data", "weights")
• func
  • オリジナルの実装がどこで定義されているか、パッケージ名と関数名を文字列で指定。
  • 例: c(pkg = "MASS", fun = "rlm")
• defaults
  • 引数の既定値を変更したい場合、引数の名前と値をリストで与える
  • 例: list()
    • 何も変更しない

parsnip::set_fit(
  model = "linear_reg",
  eng = "rlm",
  mode = "regression",
  value = list(
    interface = "formula",
    protect = c("formula", "data", "weights"),
    func = c(pkg = "MASS", fun = "rlm"),
    defaults = list()
  )
)

4. 予測方法定義

くどいですが、予測の定義にもmodel、eng、mode引数を指定します。 加えて

• type引数
  • 返り値の型を指定
  • 例: "numeric"
• value引数
  • 内部での挙動を名前付きリストで指定
    • pre: 前処理用関数またはNULL
    • post: 後処理用関数またはNULL
    • func: 予測に使う関数の名前
      • 例: c(fun = "predict")
    • args: tidymodelsからpredict関数に渡る引数のリスト
      • rlang::expr関数を用いることで、tidymodelsの内部変数を利用できる。
        • とりあえずobject = expr(object$fit)とnewdata = expr(new_data)はそういうものと思っておこう
      • 例: list(object = expr(object$fit), newdata = expr(new_data), type = "response")

parsnip::set_pred(
  model = "linear_reg",
  eng = "rlm",
  mode = "regression",
  type = "numeric",
  value = list(
    pre = NULL,
    post = NULL,
    func = c(fun = "predict"),
    args =
      list(
        object = rlang::expr(object$fit),
        newdata = rlang::expr(new_data),
        type = "response"
      )
  )
)

5. テスト

library(magrittr)
parsnip::linear_reg() %>% 
  parsnip::set_engine("rlm") %>% 
  parsnip::fit(mpg ~ ., data = mtcars)
## parsnip model object
## 
## Fit time:  6ms 
## Call:
## rlm(formula = mpg ~ ., data = data)
## Converged in 8 iterations
## 
## Coefficients:
## (Intercept)         cyl        disp          hp        drat          wt 
## 17.82250038 -0.27878615  0.01593890 -0.02536343  0.46391132 -4.14355431 
##        qsec          vs          am        gear        carb 
##  0.65307203  0.24975463  1.43412689  0.85943158 -0.01078897 
## 
## Degrees of freedom: 32 total; 21 residual
## Scale estimate: 2.15

ハマりポイント

いきなりですが、以下の点でハマりました。

1. interfaceがtidymodelsと整合しないため、オリジナルパッケージでラッパー関数を定義する必要がある
2. predictメソッドが定義されていないため、自前で用意する必要がある
3. “predictor_indicators”がどうのこうのと言ってエラーを吐くので、parsnip::set_encoding関数が必要になった
4. パッケージ読み込み時にモデルをtidymodelsに登録するには、.onLoad関数を使う
5. パッケージ開発中にdevtools::loadやdevtools::documentすると、その度に.onLoadが走って登録済みのモデルを上書きできないとtidymodelsに怒られる。parsnip::show_engines("linear_reg")を見て、登録済みなら警告しつつ、登録をスキップする処理を入れておくと良い。

実装

説明性重視で、パッケージのコードに比べ、一部簡略化してます。

horseshoe <- function(x, y, ...) {
  structure(
    horseshoe::horseshoe(y = y, X = x, ...),
    class = c("horseshoe", "list")
  )
}

predict.horseshoe <- function(object, newdata, ...) {
  newdata %*% object$BetaHat
}

テスト

\(y = 3x + 3\)に従うデータを乱数で用意し、\(y = a_0 + a_1 x + a_2 x^2 + a_3 x^3\)で回帰してみる。

library(tidyhorseshoe)
library(magrittr)
set.seed(0)

x <- runif(4, 0, 5)
y <- 3 * x + rnorm(length(x), 0, 0.1) + 3
d <- data.frame(y, x = x, x2 = x * x, x3 = x * x * x)

fitted <- parsnip::linear_reg() %>% 
  parsnip::set_engine("horseshoe") %>% 
  parsnip::fit(y ~ ., data = d)
## [1] 1000
## [1] 2000
## [1] 3000
## [1] 4000
## [1] 5000
## [1] 6000

fitted$fit$BetaMedian
## [1]  2.384374255  2.966487234  0.051657307 -0.008521586

cbind(y, predict(fitted, d))
##           y  .pred_V1
## 1 16.583438 16.584077
## 2  7.109873  6.768883
## 3  8.623323  8.403944
## 4 11.438805 11.547175

よさそう。

ちなみに、モデル式を y ~ 0 + .とすると、y切片が0になる。 parsnip::set_encoding関数でcompute_intercept = FALSEをよくわからず指定してるけれど、とりあえずはこれでよさそう。

fitted_without_intercept <- parsnip::linear_reg() %>% 
  parsnip::set_engine("horseshoe") %>% 
  parsnip::fit(y ~ 0 + ., data = d)
## [1] 1000
## [1] 2000
## [1] 3000
## [1] 4000
## [1] 5000
## [1] 6000

fitted_without_intercept$fit$BetaMedian
## [1]  4.87120134 -0.14234623 -0.01764793