The One with ...

思いついたことや作業メモなどを公開しています

探索的分析の自動化実験

内容

データフレームを突っ込むと,データに含まれる全変数の中から1つを順番に応答変数として選び,残りの変数に回するような分析を自動的にやってくれるコードを書いてみた.

「本格的な分析の前に,変数間の関係を探索的・網羅的に知りたい」場合に使うことを想定している. 「こういう頭を使わない機械的な分析はやってはいけない」という例を示すための実験でもある.

目的

データがn個の変数{X_1, X_2, \ldots, X_n}を含む場合に,

{X_1}が応答変数で,それ以外の全ての変数が説明変数
{X_2}が応答変数で,それ以外の全ての変数が説明変数
・・・
{X_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するコードが別途必要になる.ちなみに上の二つのコードはこちらからの引用である.

stackoverflow.com

上記のコードを元に,自動データ分析関数のプロトタイプを作ってみた.

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の使い方は意味が無いし,「有意性を見つけてから理屈を考える」というアドホックな分析は,研究としてやってはいけない. 上記の関数からは,理論的に無意味なアウトプットもたくさん出てくるので,使用はリコードの確認や 探索的分析にとどめよう.