Rのモデル式からの情報取得
Rを使い始めると初日からお世話になるy ~ x + zのようなモデル式ですが、意外に扱い方に習熟しづらいものです。関数にモデル式を入れるのは簡単ですが、モデル式を操作するのに慣れていない人は多いでしょう。しかし、パッケージの作成などでは必要になる知識です。慣れてしまえば簡単なので、存在を覚えておきましょう。
1 モデル式の作り方
lmなどのコマンドに渡すモデル式を作る方法はOLSの説明でざっと説明しました。?formulaにも同様のことが書いてあります。
2 モデル式からの情報取得
つくったモデル式は、lmなどのコマンド内で情報取得に使われるわけですが、何をどうやって取得するのか説明されることは少ないです。パッケージ作成などで汎用性の高い関数をつくるときには有用なので、ざっと確認していきましょう。
2.1 データセット
例に使うデータフレームdf01は以下で、100行あります。
| g | x1 | x2 | d |
|---|---|---|---|
| 0 | 0.4423369 | 1.2189649 | c |
| 2 | 0.0831099 | 0.0734796 | a |
| 0 | 0.6146112 | 1.0154364 | b |
| 2 | 0.6228064 | 0.6251571 | c |
| 0 | 0.5373650 | 0.8203078 | b |
2.2 変数名の取得
モデル式g ~ x1 + x2 + dから、被説明変数（従属変数）と説明変数（独立変数）を抜き出してみましょう。
getvars <- function(frml){
# 全ての変数名を取得
avars <- all.vars(frml)
# モデル式をリスト化
lst <- as.list(frml)
# リストは~と左辺と右辺に分解されるので、右辺の全変数が説明変数となる
evars <- all.vars(lst[[length(lst)]])
# 被説明変数を選択
dvars <- avars[!avars %in% evars]
# リストに整理して戻す
list(all = avars, dependent = dvars, explanatory = evars)
}
getvars(g ~ x1 + x2 + d)$all
[1] "g" "x1" "x2" "d"
$dependent
[1] "g"
$explanatory
[1] "x1" "x2" "d" かなり複雑なモデル式でも、変数名だけを抜き出せます。
$all
[1] "y1" "y2" "x1" "x2"
$dependent
[1] "y1" "y2"
$explanatory
[1] "x1" "x2"2.3 モデル式に必要となる変数をまとめる
model.frameを使うとモデル式で推定するのに必要な変数をデータフレームにまとめることができます。
なお、属性はダミー変数に展開されません。I(...)は実際の値に展開されますが、クロス項は省略されています。
| g | I(x1^2) | x1 | x2 | d |
|---|---|---|---|---|
| 0 | 0.195661…. | 0.4423369 | 1.2189649 | c |
| 2 | 0.006907…. | 0.0831099 | 0.0734796 | a |
| 0 | 0.377746…. | 0.6146112 | 1.0154364 | b |
| 2 | 0.387887…. | 0.6228064 | 0.6251571 | c |
| 0 | 0.288761…. | 0.5373650 | 0.8203078 | b |
属性情報が付与されているのでtermsで参照することができます。lmやglmの結果オブジェクトにmodel.frameでつくったデータフレームが入っていますが、予測値を作るときなどに使えます。
g ~ I(x1^2) + x1 * x2 + d
attr(,"variables")
list(g, I(x1^2), x1, x2, d)
attr(,"factors")
I(x1^2) x1 x2 d x1:x2
g 0 0 0 0 0
I(x1^2) 1 0 0 0 0
x1 0 1 0 0 1
x2 0 0 1 0 1
d 0 0 0 1 0
attr(,"term.labels")
[1] "I(x1^2)" "x1" "x2" "d" "x1:x2"
attr(,"order")
[1] 1 1 1 1 2
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: 0x6447ff98be90>
attr(,"predvars")
list(g, I(x1^2), x1, x2, d)
attr(,"dataClasses")
g I(x1^2) x1 x2 d
"ordered" "numeric" "numeric" "numeric" "character" ダミーが含まれるモデルで、新たにnewdataを与えて予測値をつくるときなどで、character型を元のデータフレームの列には無い要素もlevelsに含むfactor型にしたい場合は、xlevオプションが使えます。
2.4 欠損値を含める
標準ではoptions(na.action = "na.omit")になっていて、model.frameはNAを含む行を除外します。欠損値処理がlist-wise/pair-wiseになってしまうわけですが、補定する場合は都合が悪いでしょう。model.frameの引数にna.action = na.passをつけるか、options(na.action = "na.pass")とグローバルオプションを変更すれば、NAを含む行も含むデータフレームが作成されます。
2.5 モデル式から行列を作成する
同様に行列にすることもできます。factor型やcharacter型の列はダミー変数に展開され、クロス項も計算されます。一方、属性情報は付与されないのでtermsは使えません。
| (Intercept) | I(x1^2) | x1 | x2 | db | dc | x1:x2 |
|---|---|---|---|---|---|---|
| 1 | 0.1956619 | 0.4423369 | 1.2189649 | 0 | 1 | 0.5391932 |
| 1 | 0.0069072 | 0.0831099 | 0.0734796 | 0 | 0 | 0.0061069 |
| 1 | 0.3777470 | 0.6146112 | 1.0154364 | 1 | 0 | 0.6240986 |
| 1 | 0.3878878 | 0.6228064 | 0.6251571 | 0 | 1 | 0.3893518 |
| 1 | 0.2887611 | 0.5373650 | 0.8203078 | 1 | 0 | 0.4408047 |
この命名規則が仕様かは明確な記述を見つけられていませんが、factor型やcharacter型の列をダミー変数に展開するとき、列名が属性名（ここではd）とその値になりうるカテゴリー名の合成になります。カテゴリー名はcharacter型ではデータとして入っている値からつくられますが、factor型の場合はlevelsの値が参照されます。なお、切片項がある場合は、最初の変数は省略されます。
| (Intercept) | x1 | x2 | db | dc |
|---|---|---|---|---|
| 1 | 0.4423369 | 1.2189649 | 0 | 1 |
| 1 | 0.0831099 | 0.0734796 | 0 | 0 |
| 1 | 0.6146112 | 1.0154364 | 1 | 0 |
| 1 | 0.6228064 | 0.6251571 | 0 | 1 |
| 1 | 0.5373650 | 0.8203078 | 1 | 0 |
先ほどmodel.frameでxlevを設定してつくったデータフレームではダミーの列が増えています。
| (Intercept) | x1 | x2 | db | dc | dZ |
|---|---|---|---|---|---|
| 1 | 0.4423369 | 1.2189649 | 0 | 1 | 0 |
| 1 | 0.0831099 | 0.0734796 | 0 | 0 | 0 |
| 1 | 0.6146112 | 1.0154364 | 1 | 0 | 0 |
| 1 | 0.6228064 | 0.6251571 | 0 | 1 | 0 |
| 1 | 0.5373650 | 0.8203078 | 1 | 0 | 0 |