<< 目次へ
最終更新日： 2010.5.30
変数サンプリングによるモデル選択と欠損値補完
+ 個別関数（拡張子 .R で保存されたファイル） |- rslm() # RSLMを実行する |- rslm.simulate() # RSLMによって最適なモデルを探索する。rslm()が含まれる。 |- mlri() # 欠損値の補完を行う。データ行列にNAが含まれないもの用（シミュレーション用）。 |- mlri.simulate() # mlri()を繰り返し実行させる関数。 |- mlri2() # 欠損値の補完を行う。データ行列にNAが含まれているもの用（現場用）。 |- atehame() # 判別的中率を返すための関数。mlri類の関数内で必要。 + 読み込み用Rファイル（拡張子 .R で保存されたファイル） |- rslm.load.R # rslm()とrslm.simulate()が含まれている。読み込みようのファイル。 |- mlri.load.R # シミュレーション用のmlri()を実行するための全ての関数が含まれている。 |- mlri2.load.R # 現場用のmlri2()を実行するための全ての関数が含まれている。
RSLM（Random Sampling Linear Model）
変数をランダムサンプリングして、そのサンプリングされた変数を説明変数として重回帰分析を行う。これ自体はあまり意味のない作業だが、繰り返し行うことで効率よく決定係数の高いモデルを探索することができる。これを実行するためには読み込み用Rファイル rslm.load.R を作業ディレクトリなどに保存しておき、[ファイル] -> [Rファイルのソースコードを読み込み]から読み込めばよい。使用するデータセットは missing.txt であるが、これには最終4行にNAが含まれた行があるのでそれを除いたものを用意しておく必要がある。
> missing <- read.table("missing.txt", header=F) # missing.txtの読み込み
> dim(missing) # データ行列の次元を確認する
[1] 104 12
> missing[1:2,] # 最初の2行目まで表示
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
1 1 2 3 2 2 4 1 1 1 1 1 1
2 2 2 2 3 2 1 2 1 1 1 2 2
> rslm(dat) # 第1引数には欠損値を除いたデータ行列を指定する
$MODEL モデル式
V12 ~ V6 + V11 + V9
$RESULT 解析結果
Call:
lm(formula = model, data = dat)
Coefficients:
(Intercept) V6 V11 V9
-0.01663 0.15709 0.26769 0.19213
$ARS 自由度調整済み決定係数
[1] 0.4627193
rslm()の結果はリスト形式で返される。それにはMODEL（モデル式）、RESULT（lm()の解析結果）、ARS（自由度調整済み決定係数）が含まれており、rslm()の返す結果をオブジェクト（変数）に代入しておけば、自由に取り出すことができる。なお、返される結果はrslm()を実行するたびに異なるので注意した方がよい。
> result <- rslm(dat) > names(result) [1] "MODEL" "RESULT" "ARS" > result$MODEL V12 ~ V5 + V2 > result$ARS [1] 0.2363030
最初に述べたように、rslm()を繰り返し実行し、モデル更新のたびに決定係数の高いモデルを選んでいけば、最終的には繰り返し数の中で最も高い決定係数を有していたモデルを得ることができる。これを行うための関数がrslm.simulate()である。第1引数には最終列に応答変数をおいたデータ行列を、第2引数には繰り返し数を指定する。繰り返し数は100くらいで十分である。
> rslm.simulate(dat, 100) # 100回の繰り返しでrslm()を実行させる
model.formula:
V12 ~ V4 + V10 + V8 + V11 + V3 + V5 + V2 + V6 + V1 + V9
Adj.R.Squared:
[1] 0.5260338
model.formula: 選ばれたモデル
V12 ~ V9 + V4 + V6 + V1 + V11 + V5
Adj.R.Squared: 自由度調整済み決定係数
[1] 0.5328212
Complete!! 完了の合図
Call:
lm(formula = model, data = dat)
Coefficients:
(Intercept) V9 V4 V6 V1 V11
-0.36698 0.18061 -0.13496 0.06069 0.27606 0.20505
V5
0.17777
RSLMと多項ロジットモデルによる欠損値補完
応答変数に欠損値がある場合、あらかじめRSLMによって決定係数の高いモデルを探し、そこから得られる予測値を説明変数として多項ロジットモデルを当てはめる。これによって得られる予測確率に基づいて、適当な補完値を割り当てるという方法を適用してみる。現場で使えるもの、つまり実際に欠損値のあるデータセットに対して使えるmlri2()と、シミュレーション用に値を伏せて、それをどれだけリカバーできるかを試すために使えるmlri()の2つがある。それぞれ mlri2.load.R と mlri.load.R を読み込めば必要な関数が全てR側にロードされる。
なおどちらにも判別的中率を返す関数 atehame() が含まれているが、mlri2()を実行した際には欠損値に対する判別的中率は表示されない（シミュレーション用と違って、値を伏せているわけではないのだから確認の使用がないというわけ）。それから mlri.load.R には mlri.simulate() が含まれているが、これは判別的中率が繰り返し数や欠損値数によってどのように変化していくかをシミュレートするためのものである。
まず欠損値が含まれたデータセットmissingを使って、mlri2()を実行した結果を以下に示す。ただし、出力が膨大なものになるので、sink("out.txt")などとしてテキストファイルに結果を出力した方が良いかもしれない。sink()は出力先を変更するための関数だが、引数なしでsink()とコマンドすればR側に戻ってくる。
それからmlri()もmlri2()も関数内でVGAMパッケージに含まれるvglm()を用いている。VGAMが呼び出されていない場合は自動的に呼び出されるようになっているが、インストールされていないとエラーが出て使えないので、あらかじめダウンロードしてインストールしておく必要がある。
第1引数には欠損値を含んだ、かつ最終列に応答変数をおいたデータ行列を、第2引数にはrslm()を繰り返す回数を指定する。この繰り返し数は100程度で十分である。ただし、rslm()の出力はRのウィンドウ（画面）上には表示されず、別途rslm.out.txtというテキストファイルとして作業ディレクトリに出力される。
> mlri2(missing, 100)
rslm.out.txtとしてワークディレクトリに出力されました
$MULTIPLE 重回帰の結果
Call:
lm(formula = model, data = dat)
Coefficients:
(Intercept) Var11 Var10 Var1 Var5 Var9
-0.35139 0.21002 0.06955 0.32321 0.20670 0.15011
Var2 Var4
-0.09936 -0.11737
$LOGIT 多項ロジットモデルの結果
Call:
vglm(formula = as.factor(ad.Y) ~ ad.lm.prd, family = multinomial)
Coefficients:
(Intercept):1 (Intercept):2 (Intercept):3 ad.lm.prd:1 ad.lm.prd:2
21.445056 17.326485 6.331958 -8.636167 -6.082249
ad.lm.prd:3
-1.633009
Degrees of Freedom: 300 Total; 294 Residual
Residual Deviance: 158.9578
Log-likelihood: -79.47891
$ANA.MR.PREDICT 欠損値を除いた部分のデータに対する重回帰の予測値
1 2 3 4 5 6 7
0.3814467 0.7973112 2.9768827 1.0668650 1.8934890 1.4058417 1.6825976
・・・ （途中省略） ・・・
99 100
2.7357154 3.0497543
$ANA.LR.PREDICT 欠損値を除いた部分のデータに対するロジットモデルの予測確率
1 2 3 4
1 0.9586833391 0.04131286 3.786910e-06 1.255668e-08
2 0.8891133868 0.11082157 6.462352e-05 4.225854e-07
3 0.0024153032 0.07872558 7.471805e-01 1.716786e-01
・・・ （途中省略） ・・・
99 0.0118198058 0.20809579 6.754142e-01 1.046702e-01
100 0.0014514518 0.05698658 7.479812e-01 1.935808e-01
$MIS.MR.PREDICT 欠損値の部分のデータに対する重回帰の予測値
101 102 103 104
1.306185 2.210559 2.999004 1.632327
$MIS.LR.PREDICT 欠損値の部分のデータに対するロジットモデルの予測確率
1 2 3 4
1 0.685033860 0.31318244 0.001757316 2.637972e-05
2 0.139979120 0.64452195 0.202206521 1.329241e-02
3 0.002071587 0.07144695 0.748236166 1.782453e-01
4 0.481378227 0.50619008 0.012121742 3.099472e-04
$ANA.RESULT 欠損値を除いた部分のデータへ予測確率に基づいて当てはめられた値
$ANA.RESULT$FITTED.VALUE
[1] 1 1 3 1 2 1 2 3 2 2 2 1 3 1 3 3 2 2 2 1 2 1 2 1 3 1 3 2 2 2 1 2 3 2 2 3 3
[38] 2 2 2 2 3 1 3 3 2 3 3 3 2 2 2 1 1 3 1 3 3 3 2 2 2 3 3 2 1 3 3 2 2 1 3 2 2
[75] 2 1 2 3 2 1 1 3 2 2 2 3 2 2 2 2 3 2 2 2 2 3 3 2 3 3
$ANA.RESULT$TABLE 実測値と当てはめ値との差を集計したもの
0 1
63 37
$ANA.RESULT$PROBABILITY 上記集計値の割合
0 1
0.63 0.37
$MIS.RESULT 欠損値へ補完されるべき値（最終的な補完値）
$MIS.RESULT$FITTED.VALUE
[1] 1 2 3 2
$MIS.RESULT$TABLE
character(0)
$MIS.RESULT$PROBABILITY
numeric(0)
続いて欠損値のないデータに対して、任意の応答変数の値を伏せて（欠損値とみなして）、その値をどれほどリカバーできているかmlri()を用いて実行してみる。得られる出力は上記のmlri()と同じであるが、引数が少し異なる。第1引数にはデータ行列、第2引数には伏せる値の数、第3引数にrslm()の繰り返し数を指定する。
mlri()では伏せた値に対して補完値がどれだけ一致していたかという判別的中率が出力されていることに注意（先ほどは出力されていなかった）。
> mlri(dat, 10, 100)
rslm.out.txtとしてワークディレクトリに出力されました
$MULTIPLE
Call:
lm(formula = model, data = dat)
Coefficients:
(Intercept) Var7 Var5 Var9 Var6 Var1
-0.44665 0.15197 0.17524 0.13742 0.06752 0.25866
Var4 Var11
-0.12549 0.18156
$LOGIT
Call:
vglm(formula = as.factor(ad.Y) ~ ad.lm.prd, family = multinomial)
Coefficients:
(Intercept):1 (Intercept):2 (Intercept):3 ad.lm.prd:1 ad.lm.prd:2
21.698288 17.173983 7.367337 -8.802278 -5.987735
ad.lm.prd:3
-2.010782
Degrees of Freedom: 270 Total; 264 Residual
Residual Deviance: 144.2028
Log-likelihood: -72.10139
$ANA.MR.PREDICT
1 2 3 4 5 7 9
0.6525446 0.9166816 2.9495950 0.9974802 1.6235200 1.6854375 2.1492357
94 96 97 98 99 100
1.6380696 2.5259918 2.7190041 2.3285310 2.5379219 3.2598618
$ANA.LR.PREDICT
1 2 3 4
1 0.9362522738 0.06370060 4.701392e-05 1.102760e-07
2 0.8745881392 0.12514674 2.640655e-04 1.053488e-06
3 0.0024091130 0.10528961 7.208796e-01 1.714217e-01
89 0.0286621625 0.39321724 5.236977e-01 5.442287e-02
90 0.0002732991 0.02860344 6.726338e-01 2.984894e-01
$MIS.MR.PREDICT
54 84 6 95 47 26 63
1.0541029 2.2610100 1.7474877 2.4232747 2.9937373 0.8643147 2.9408769
11 8 16
1.7901653 2.3184300 3.3273828
$MIS.LR.PREDICT
1 2 3 4
1 0.8254663397 0.17389656 0.0006337718 3.333166e-06
2 0.0976091406 0.61423358 0.2719619562 1.619533e-02
3 0.3891044837 0.57704559 0.0331470440 7.028812e-04
4 0.0499630372 0.49640319 0.4190517588 3.458201e-02
5 0.0017880636 0.08848469 0.7220819190 1.876453e-01
6 0.8899212294 0.10988981 0.0001882796 6.760683e-07
7 0.0025538070 0.10890804 0.7202438465 1.682943e-01
8 0.3585841184 0.59965516 0.0408176302 9.430920e-04
9 0.0782063596 0.57845800 0.3218254155 2.151023e-02
10 0.0001666810 0.02109598 0.6489038295 3.298335e-01
$ANA.RESULT
$ANA.RESULT$FITTED.VALUE
[1] 1 1 3 1 2 2 2 2 1 2 2 3 1 2 1 2 1 1 2 1 3 3 2 2 2 1 2 3 2 2 3 3 1 1 2 2 3
[38] 1 3 3 2 3 3 2 2 2 2 3 1 3 2 3 2 2 1 3 2 1 3 3 2 2 1 3 2 2 2 1 2 3 2 1 1 3
[75] 2 2 3 3 2 2 2 3 2 2 2 3 3 2 3 3
$ANA.RESULT$TABLE
0 1 2
56 33 1
$ANA.RESULT$PROBABILITY
0 1 2
0.62222222 0.36666667 0.01111111
$MIS.RESULT
$MIS.RESULT$FITTED.VALUE
[1] 1 2 2 2 3 1 3 2 2 3
$MIS.RESULT$TABLE
0 1
6 4
$MIS.RESULT$PROBABILITY 判別的中率
0 1
0.6 0.4
最後にmrli.simulate()の使い方を以下に示す。これの引数は第1引数にデータ行列、第2引数に伏せる値の数、第3引数にrslm()の繰り返し数、第4引数にmrli()の繰り返し数を指定する。
> mlri.simulate(dat, 10, 100, 3)
rslm.out.txtとしてワークディレクトリに出力されました
--- 1 回目の繰り返し---
0 1
0.6333333 0.3666667
0 1
0.4 0.6
rslm.out.txtとしてワークディレクトリに出力されました
--- 2 回目の繰り返し---
0 1
0.6111111 0.3888889
0 1
0.4 0.6
rslm.out.txtとしてワークディレクトリに出力されました
--- 3 回目の繰り返し---
0 1
0.6 0.4
0 1
0.5 0.5
