概要

formula オブジェクトは変数変換や交互作用項など, 多彩な表現ができる.
xgboost や glmnet では data.matrix() を併用することで formula を利用できる.

統計モデリング/機械学習で予測モデルを構築するとき, 予測性能の向上のため, しばしば変数を入れ替えたり, 変換したりといった推敲が必要となる. R ではこういうときに formula オブジェクトを使うと, いちいちデータフレームに変換後の数値を追加したり書き換えたりする必要がなくなる. glmnet や xgboost など, formula が直接使えないものでも model.matrix() 等を併用すれば可能である*1.

formula オブジェクトを解説した記事を探すと, かなり前から存在する. 例えば以下の記事.
m884.hateblo.jp

なお, 上記はタイトルが
「formula とは? (1)」となっているが, 続編は確認できなかった.

もうひとつ有用な日本語のページがある.
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/71.html

今回は, これらを記事を踏まえて formula のもう少し応用的な使い方を紹介する.

ちなみに, 以前に glmnet 限定で似たようなことを ill-identified.hatenablog.com で書いたのだが, 当時はまだ model.matrix() の挙動をよく理解していなかったのであまりきれいな書き方とは言えない. そのうちこの記事も修正しようと思う.

formula オブジェクト

formula 構文はモデルの構造を記述に便利だが,

y~x1+x2

というように最も簡単な回帰式にしか使っていない人も多いと思う. しかし, この formula にはいろいろな構文がある. 例えば,

y ~ .

ならば入力データの全変数を説明変数とする.

y ~ . - x1

ならば, 全変数から x1 を除いたものを説明変数とする.

y ~ x1 * x2

ならば, x1 と x2 と, 2変数の交互作用項が追加される*2 *3. また, log(), sqrt(), sin() などの関数も利用できる. つまり,

log(y) ~ x1 + x2

ならば片対数モデルになる. log(x+1) ~ x1 + x2 も可能. さらに, I() または poly() で多項式回帰をすることも可能である. lag() を使えば自己回帰項を追加できる. また, 従来は lm() のオプション引数で指定することが多いオフセット項やウェイト変数も, offset(), weights() を使うことで formula の中で記述できる.

これらの演算子は, () で囲むことで優先順位を変えることもできるため, かなり複雑な変数の組み合わせも可能である*4. 例えば, 以上の構文の組み合わせ技として,

log(y) ~(. + poly(x1, 3)) * z1 - z1:x1

という式を書いた場合, 従属変数は y の対数となる. 説明変数は, データフレームの全変数と x1 の3次多項式に加え, これらと z1 の交互作用項を加え, そこから z1 と x1 の交互作用項だけ除いたものになる.

以上から, formula をうまく活用できれば, 元のデータフレームに変数 x1 を2乗しただけ, 対数変換しただけ, といった列を大量に作成する必要がなくなる.

ここでは動作確認用サンプルとして, 以下のようなデータフレームを作成する.

set.seed(42)
n <- 30
df.train <- data.frame(y=rnorm(n)^2,
                       x1=rnorm(n),
                       x2=rnorm(n),
                       z1=as.factor(sample(c("A", "B", "C"),
                                           size = n, replace = T,
                                           prob = c(.4, .4, .2))
                                    ),
                       z2=as.factor(sample(c("A, "B", "C";),
                                           size = n, replace = T,
                                           prob = c(.2, .2, .6))
                                    )
                 )
df.test <- df.train[floor(n/2):n, ]
df.test$y <- NULL
df.train <- df.train[1:(floor(n/2)-1), ]

このデータフレームは x1, x2 は numeric で, z1, z2 は 3種類の値をとる factor である. 予測モデルを作るという前提なので, 訓練データ df.train とテストデータ df.test を用意している. df.test には従属変数 y は存在しない. 以降はこの2つのデータフレームがあるという前提で話をすすめる.

formula オブジェクトは線形回帰モデルを推定する lm() や, 一般化線形回帰モデルを推定する glm() など, 多くの関数で使用できる. まずは片対数モデル $\ln y=\alpha+\beta_{1}x_{1}+\beta_{2}x_{2}$ の回帰をしてみる. わざわざ新しくyの対数を作成する必要はなく, 先ほど紹介したように log(y) を使って,

(m1 <- lm(formula=log(y)~x1+x2, data=df.train))

とする. さらにテストデータにモデルの予測値を与えたい場合は, 通常通り predict() を使い,

predict(m1, df.train)

とする.

あるいは $y=\alpha+\beta_{1}x_{1}+\beta_{2}x_{1}^{2}+\beta_{3}x_{1}^{3}$ のような多項式回帰なら,

(m2 <- lm(formula=y~poly(x1, 3), data=df.train))

となる. これもやはり predict() が利用可能である. これらは lm オブジェクトが formula オブジェクトの情報を保持しているので可能なことである.

formula が使えない場合

しかし, 一部のパッケージでは formula オブジェクトに対応していないものがある. 具体的には glmnet と xgboost である. これらは formula 関数ではなく, matrix オブジェクトしか受け付けない*5. このような場合は model.matrix() を使う. この関数はデータフレームを formula の式に沿って matrix に変換してくれる. ところで, かりに特殊な変数変換を行わない場合でも, factor 型や logical 型変数を含むデータを使う場合は注意が必要である. 計算をおこなううえで, factor/logical 型はゼロ-1のバイナリ変数をカテゴリの数だけ作成することが必要になるからである. 以下イメージ.

データフレームが

z1 (factor)	flag (logical)
A	TRUE
B	TRUE
C	FALSE

となっていた場合, 以下のように変換することになる.

z1A	z1B	z1C	flagTRUE	flagFALSE
1	0	0	1	0
0	1	0	1	0
0	0	1	0	1

この変換後の数値化された行列はデザイン行列と呼ばれる. lm()等の関数は処理の中にこのデザイン行列への変換処理が組み込まれているので書く必要がなかったのだ.

model.matrix() はこのようなデザイン行列への変換処理もおこなってくれる. たとえば,

model.matrix(log(y)~x1+x2+z1*z2, df.train)

のようにすれば, 式の右辺, つまり $x_{1}$ , $x_{2}$ に加え, $z_{1}$ , $z_{2}$ とその交差項のカテゴリ別のバイナリ変数を含むデザイン行列ができる.

model.matrix() を使えばこのようなデザイン行列を作成するのが楽になるが, 従属変数の変換はしてくれない. どうせならば formula を与えればすべて自動でやってくれるようにすれば楽なので, xgboost を例にして, 以下のように formula とデータフレームから xgb.Dmatrix オブジェクトを作成する関数を定義する.

require(xgboost) ver. 0.6-4
# 片対数モデルを定義
model <- log(y)~x1+x2+z1*z2

# 変換用の関数定義
xgb.DMatrix2 <- function(formula, data){
  xgb.DMatrix(data = model.matrix(formula, data),
              label = model.response(model.frame(formula, data)))
}

# xgboost で訓練
train <- xgb.train(data = xgb.DMatrix2(model, df.train),
                   params = list(objective=&#39;reg:linear&#39;),
                   nrounds=10)

# 訓練データに対する予測
predict(train, newdata=xgb.DMatrix2(model, df.train))

ここでは, model.response() を使って formula の左辺で評価した値をラベル値に与えている. この model.response() は model.frame オブジェクトのみを引数に取るため, model.frame() 関数を挟んでいる.

さらにテストデータに対して予測値を与える場合は, テストデータのデザイン行列が必要になるが, テストデータには従属変数がふくまれないためエラーとなる. この場合は以下のようにしてデザイン行列を作成できる.

model.matrix(model[-2], df.test)

別解として,

model.matrix(delete.response(terms(model)), df.test)

とすることもできる. この話から分かるように, 先に定義した xgb.DMatrix2() にテストデータを与えるとエラーが発生するので, より便利にするため以下のように改良する.

# 変換用の関数定義
xgb.DMatrix2 <- function(formula, data){
  if(prod(all.vars(model[-3]) %in% colnames(data)) == 1){
    # data contains dep. variables
    return(xgb.DMatrix(data = model.matrix(formula, data),
                       label = model.response(model.frame(formula, data)))
           )
  }
  else{
    # data does not contain dep. variables
    return(xgb.DMatrix(data=model.matrix(formula[-2], data)))
  }
}

# テストデータで予測
predict(train, newdata=xgb.DMatrix2(model, df.test))

これで, 従属変数を含むデータを与えた場合は訓練データとみなし, そうでない場合はテストデータとみなして処理を場合分けするようになった. ただし, このやり方には2点, 注意すべきことがある. 1つは, テストデータに y という名前の変数があるものの, 値がすべて NA であるようなデータを与えてしまった場合である. このとき, 欠損値が含まれている行は除外されてしまうので, 上記の関数では0行のマトリクスを返してしまうから, y そのものを消す必要がある. もう1つは, この方法では予測値を従属変数の変換後の値で返してしまう, ということである. 今回の例では model に片対数モデルが定義されているから, predict() で返されるのは y の予測値ではなく対数の予測値であり, 正しくは exp(predict(...)) と書く必要がある.

*1:ただし, caret の train() を併用するという手もある. 少し前のバージョンではうまくいかなかった記憶があるのだが, 現時点の最新版 ver. 6.0-76 ではうまくいくようである. もう少し早くにこの記事を書いておけばよかった.

*2:相互作用項だけがほしいならコロン : を代わりに使う.

*3:2変数がいずれも factor でないなら, 単なる積値になる.

*4:蛇足だが, SAS の reg プロシジャ等でもこれに似た構文が使える.

*5:先に書いたように caret パッケージと併用すれば formula を利用することができる.