読者です 読者をやめる 読者になる 読者になる

昨日もらったコメントを元に昔書いたのをsapplyで書きなおしてみるテスト

R

昔書いた全ての組み合わせを作り出す的なスクリプトをsapply関数を使ってスマートに書いてみた。最初の部分は同じ。

car <- c("height","weight","length","eu")
grid <- eval(parse(text=paste("expand.grid(",paste(car,collapse=" = c(0,1),")," = c(0,1))")))
#追記しました
#上のほうが、変更しなくてよいので、楽チン
#grid <- expand.grid(height = c(0,1), weight = c(0,1), length = c(0,1),eu = c(0,1))

reg <- c()

for(i in 2:length(grid[,1])){
  reg[i-1] <- list(car[grid[i,]==1])
}

で、昨日のid:arboresさんにコメントをいただいたところを活用させてもらって、sapply関数でこんなことをやってみる。コメントにもあるようにpaste関数はベクトルを引数に取ることができるのを利用している*1

> sapply(reg,function(x){paste(x,collapse="+")})
 [1] "height"                  "weight"                 
 [3] "height+weight"           "length"                 
 [5] "height+length"           "weight+length"          
 [7] "height+weight+length"    "eu"                     
 [9] "height+eu"               "weight+eu"              
[11] "height+weight+eu"        "length+eu"              
[13] "height+length+eu"        "weight+length+eu"       
[15] "height+weight+length+eu"

でもって、このベクトルにさらにpasteをかますこともできるので、こんなことができる。

> paste("price~",sapply(reg,function(x){paste(x,collapse="+")}),sep="")
 [1] "price~height"                  "price~weight"                 
 [3] "price~height+weight"           "price~length"                 
 [5] "price~height+length"           "price~weight+length"          
 [7] "price~height+weight+length"    "price~eu"                     
 [9] "price~height+eu"               "price~weight+eu"              
[11] "price~height+weight+eu"        "price~length+eu"              
[13] "price~height+length+eu"        "price~weight+length+eu"       
[15] "price~height+weight+length+eu"

あとは昔やったように上のやつをevalして、lmとかをすればよい。例えばこんなの。Mapを使ってやるぜ、ふひひ。16個のモデルとか多かったので、tailで最後の二つだけを出力してみた。

> Map(function(x){eval(parse(text=x))},tail(paste("summary(lm(price~",sapply(reg,function(x){paste(x,collapse="+")}),"))",sep=""),n=2))
$`summary(lm(price~weight+length+eu))`

Call:
lm(formula = price ~ weight + length + eu)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.31033 -0.13538 -0.05015  0.16090  0.47192 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.03886    0.55786   0.070    0.947
weight       0.03000    0.04382   0.685    0.519
length       0.02124    0.01147   1.852    0.113
eu          -0.01007    0.25529  -0.039    0.970

Residual standard error: 0.3039 on 6 degrees of freedom
Multiple R-squared: 0.3831,	Adjusted R-squared: 0.0746 
F-statistic: 1.242 on 3 and 6 DF,  p-value: 0.3743 


$`summary(lm(price~height+weight+length+eu))`

Call:
lm(formula = price ~ height + weight + length + eu)

Residuals:
       1        2        3        4        5        6        7        8 
-0.07543 -0.29616  0.06217  0.01760  0.29182  0.12774  0.19289  0.05661 
       9       10 
-0.09804 -0.27921 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.274561   0.498632  -0.551   0.6056  
height      -0.150046   0.080506  -1.864   0.1214  
weight       0.091324   0.049417   1.848   0.1239  
length       0.025314   0.009897   2.558   0.0508 .
eu           0.503984   0.349600   1.442   0.2090  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.2557 on 5 degrees of freedom
Multiple R-squared: 0.636,	Adjusted R-squared: 0.3448 
F-statistic: 2.184 on 4 and 5 DF,  p-value: 0.2071 

for文を取り除くことができて、大変スマートです。よかったですね。

疑問

上のmap関数を使ったものと下のsapplyは同じことをしているように見えるのだけど、出力される結果が違う。どういうことなんだろう?

sapply(paste("lm(price~",sapply(reg,function(x){paste(x,collapse="+")}),")",sep=""),function(x){eval(parse(text=x))})

出力される結果としては次のようなもの。

              lm(price~weight+length+eu) lm(price~height+weight+length+eu)
coefficients  Numeric,4                  Numeric,5                        
residuals     Numeric,10                 Numeric,10                       
effects       Numeric,10                 Numeric,10                       
rank          4                          5                                
fitted.values Numeric,10                 Numeric,10                       
assign        Integer,4                  Integer,5                        
qr            List,5                     List,5                           
df.residual   6                          5                                
xlevels       List,0                     List,0                           
call          Expression                 Expression                       
terms         Expression                 Expression                       
model         List,4                     List,5                           

*1:偉そうに書いてますが、コメントを頂くまで知らなかった。id:arbores++