Rのモデル式からの情報取得

Rを使い始めると初日からお世話になるy ~ x + zのようなモデル式ですが、意外に扱い方に習熟しづらいものです。関数にモデル式を入れるのは簡単ですが、モデル式を操作するのに慣れていない人は多いでしょう。しかし、パッケージの作成などでは必要になる知識です。慣れてしまえば簡単なので、存在を覚えておきましょう。

1 モデル式の作り方
2 モデル式からの情報取得

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"

かなり複雑なモデル式でも、変数名だけを抜き出せます。

getvars(y1 + y2 ~ I(x1^2) + x1*x2 + I(x2^2))

$all
[1] "y1" "y2" "x1" "x2"

$dependent
[1] "y1" "y2"

$explanatory
[1] "x1" "x2"

2.3 モデル式に必要となる変数をまとめる

model.frameを使うとモデル式で推定するのに必要な変数をデータフレームにまとめることができます。なお、属性はダミー変数に展開されません。I(...)は実際の値に展開されますが、クロス項は省略されています。

mf <- model.frame(g ~ I(x1^2) + x1*x2 + d, data = df01)

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でつくったデータフレームが入っていますが、予測値を作るときなどに使えます。

terms(mf)

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: 0x000000002d3f3b20>
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オプションが使えます。

mf <- model.frame(g ~ x1 + x2 + d, data = df01, xlev = list(d = c("a", "b", "c", "Z")))

2.4 モデル式から行列を作成する

同様に行列にすることもできます。factor型やcharacter型の列はダミー変数に展開され、クロス項も計算されます。一方、属性情報は付与されないのでtermsは使えません。

mm <- model.matrix(g ~ I(x1^2) + x1*x2 + d, data = df01)

(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の値が参照されます。なお、切片項がある場合は、最初の変数は省略されます。

X <- model.matrix(terms(mf), df01)

(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を設定してつくったデータフレームではダミーの列が増えています。

X <- model.matrix(terms(mf), mf)

(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

R

R設定

Rプログラミング

Rで統計解析

Rグラフィックス

R拡張