探索的分析の自動化実験
内容
データフレームを突っ込むと,データに含まれる全変数の中から1つを順番に応答変数として選び,残りの変数に回するような分析を自動的にやってくれるコードを書いてみた.
「本格的な分析の前に,変数間の関係を探索的・網羅的に知りたい」場合に使うことを想定している. 「こういう頭を使わない機械的な分析はやってはいけない」という例を示すための実験でもある.
目的
データがn個の変数を含む場合に,
が応答変数で,それ以外の全ての変数が説明変数
が応答変数で,それ以外の全ての変数が説明変数
・・・
が応答変数で,それ以外の全ての変数が説明変数
という分析を自動化したい.
コードの作成
基本的なアイデアは次の2つのステップに区分できる.
一つは,テキストで式を自動生成してから,fomula
オブジェクトに変換してlm
に突っ込むという方法である.
これは,R
の型キャストを利用した方法であり,オブジェクトの使い方がやや邪道な気もするが,型を意識した
コードの理解には適している.
xnam <- paste("x", 1:25, sep="") fmla <- as.formula(paste("y ~ ", paste(xnam, collapse= "+")))
もう一つは,応答変数を先頭列にもってきて,データフレームごと,lm
に突っ込むというワイルド(?)なやり方である.
d2 <- data.frame(y,x1,x2,x3) lm(d2) Call: lm(formula = d2)
この場合,データフレームの先頭列をrotateするコードが別途必要になる.ちなみに上の二つのコードはこちらからの引用である.
上記のコードを元に,自動データ分析関数のプロトタイプを作ってみた.
analyzer <- function(dataframe){ n <- length(dataframe)#data.frame内の変数の個数を取得 no.col <- 1:n #ローテーション用ベクトル n1 <- n-1 #ループエンド用 for(i in 0:n1){ x <- no.col[1] # print(x) no.col[1] <- no.col[i+1] no.col[i+1] <- x #print(no.col) dataframe0 <- dataframe[,no.col] out = lm(formula = dataframe0) print(summary(out)) }}
計算例
なにかテキトーなデータ行列を使って試してみる.
dat0 = data.frame(swiss)#data行列の定義 head(dat0) #中身確認 analyzer(dat0)
アウトプットはこんな感じになる.
Call: lm(formula = dataframe0) Min 1Q Median 3Q Max -15.2743 -5.2617 0.5032 4.1198 15.3213 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 66.91518 10.70604 6.250 1.91e-07 *** Agriculture -0.17211 0.07030 -2.448 0.01873 * Examination -0.25801 0.25388 -1.016 0.31546 Education -0.87094 0.18303 -4.758 2.43e-05 *** Catholic 0.10412 0.03526 2.953 0.00519 ** Infant.Mortality 1.07705 0.38172 2.822 0.00734 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 7.165 on 41 degrees of freedom Multiple R-squared: 0.7067, Adjusted R-squared: 0.671 F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10 Call: lm(formula = dataframe0) Min 1Q Median 3Q Max -29.958 -9.107 1.287 10.805 24.504 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 135.1861 22.7578 5.940 5.26e-07 *** Fertility -0.7410 0.3027 -2.448 0.018727 * Examination -0.9627 0.5117 -1.881 0.067050 . Education -1.4908 0.4119 -3.619 0.000804 *** Catholic 0.1674 0.0762 2.197 0.033745 * Infant.Mortality -0.3609 0.8637 -0.418 0.678268 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 14.87 on 41 degrees of freedom Multiple R-squared: 0.618, Adjusted R-squared: 0.5715 F-statistic: 13.27 on 5 and 41 DF, p-value: 1.045e-07 Call: lm(formula = dataframe0) Residuals: Min 1Q Median 3Q Max -9.0429 -2.8491 -0.5308 3.0459 6.8874 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 24.57200 8.23969 2.982 0.00480 ** Fertility -0.09523 0.09371 -1.016 0.31546 Agriculture -0.08254 0.04387 -1.881 0.06705 . Education 0.33037 0.12858 2.569 0.01392 * Catholic -0.06778 0.02108 -3.215 0.00254 ** Infant.Mortality 0.09730 0.25297 0.385 0.70250 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.353 on 41 degrees of freedom Multiple R-squared: 0.7346, Adjusted R-squared: 0.7022 F-statistic: 22.7 on 5 and 41 DF, p-value: 7.624e-11
中略.
このように,関数を実行するとlmの結果がずらずら出てくる.
さて,以上の結果をさらにstep
関数に渡してモデル選択まで自動でやろうとすると,エラーが出てしまう.
どうやら
formula = dataframe0
の部分が,step関数のオブジェクトとして不適切らしい.
しかたないので,一旦データフレームの列をrotateしてから,順番の入れ替わった変数名リストを取り出し,
テキスト型でy ~ x1 + x2 + ... + xn
の形にしてからformula
にキャストする,という荒技(?)
を使ってみる.ついでにlmの実行の度に応答変数が変わるので,名前が分かるようにprintしておく.
#step関数を使い, AIC基準でモデル選択まで実行 analyzer101 <- function(dataframe){ n <- length(dataframe)#data.frame内の変数の個数を取得 no.col <- 1:n#ローテーション用ベクトル n1 <- n-1#ループエンド用 for ( i in 0:n1 ){ x <- no.col[1] # print(x) no.col[1] <- no.col[i+1] no.col[i+1] <- x #ベクトルのrotate #print(no.col) dataframe0 <- dataframe[,no.col]#データフレーム列のrotate variable = names(dataframe0) xnam = variable[2:n] #説明変数のラベル ynam = variable[1] #oucomeのラベル model = as.formula(paste(paste(ynam,"~") , paste(xnam, collapse= "+"))) out1 = lm(model,data=dataframe0) print(paste("Outcome is", names(dataframe0)[1])) out2 <- step(out1,trace=1)#AIC基準でモデル選択 print(out2) }}
結果を確認してみよう. データフレームはさっきと同じ.
> analyzer101(dat0) [1] "Outcome is Fertility" Start: AIC=190.69 Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC - Examination 1 53.03 2158.1 189.86 <none> 2105.0 190.69 - Agriculture 1 307.72 2412.8 195.10 - Infant.Mortality 1 408.75 2513.8 197.03 - Catholic 1 447.71 2552.8 197.75 - Education 1 1162.56 3267.6 209.36 Step: AIC=189.86 Fertility ~ Agriculture + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC <none> 2158.1 189.86 - Agriculture 1 264.18 2422.2 193.29 - Infant.Mortality 1 409.81 2567.9 196.03 - Catholic 1 956.57 3114.6 205.10 - Education 1 2249.97 4408.0 221.43 Call: lm(formula = Fertility ~ Agriculture + Education + Catholic + Infant.Mortality, data = dataframe0) Coefficients: (Intercept) Agriculture Education Catholic Infant.Mortality 62.1013 -0.1546 -0.9803 0.1247 1.0784 [1] "Outcome is Agriculture" Start: AIC=259.3 Agriculture ~ Fertility + Examination + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC - Infant.Mortality 1 38.59 9101.4 257.50 <none> 9062.9 259.30 - Examination 1 782.31 9845.2 261.20 - Catholic 1 1066.74 10129.6 262.53 - Fertility 1 1324.81 10387.7 263.72 - Education 1 2895.46 11958.3 270.33 Step: AIC=257.5 Agriculture ~ Fertility + Examination + Education + Catholic Df Sum of Sq RSS AIC <none> 9101.4 257.50 - Examination 1 816.46 9917.9 259.54 - Catholic 1 1076.03 10177.5 260.76 - Fertility 1 1947.58 11049.0 264.62 - Education 1 3091.51 12192.9 269.25 Call: lm(formula = Agriculture ~ Fertility + Examination + Education + Catholic, data = dataframe0) Coefficients: (Intercept) Fertility Examination Education Catholic 132.6210 -0.7989 -0.9802 -1.5193 0.1681
以下,略.
いまのところnumericなデータにしか使えないが,気が向いたらfactorが応答変数になった場合は自動的に 多項ロジットを実行する機能などを追加してみる.
繰り返すが,こんなAICの使い方は意味が無いし,「有意性を見つけてから理屈を考える」というアドホックな分析は,研究としてやってはいけない. 上記の関数からは,理論的に無意味なアウトプットもたくさん出てくるので,使用はリコードの確認や 探索的分析にとどめよう.