The One with ...

研究用のメモなどを公開しています

データフレームを突っ込むと,データに含まれる全変数を応答変数にした分析を自動的にやってくれるコードを書いてみた.

具体的に言えば,データが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が応答変数になった場合は自動的に 多項ロジットを実行する機能などを追加してみる.

なお上記の関数からは,理論的に無意味なアウトプットもたくさん出てくるので,使用はリコードの確認や 探索的分析にとどめましょう.