StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

Pythonのリスト内包表記みたいなのをRで書く方法(とおまけ)

Rの小技

 Rでデータ分析するときに最近よくやる書き方があるのですが、よくやるといいながら1か月とか間が開くと忘れているので、メモしておきます。
 2個あるのですが、1個目は将来ネットで検索して役に立ててくれる人がいるかもしれないので、1個目の小技をこのエントリのタイトルにして2個目をおまけ扱いにしました。


 私はべつに、Rにもプログラミングにも詳しい人間ではないので、妥当なやり方かどうかは知りません。
 というか、もっと良いやり方あればご教示いただけると・・・。(とくに2個目)
 
 

1.Pythonのリスト内包表記みたいなやつをRで書く

 Pythonだと、リストに入っている要素について、連続的に何かの処理をして、その結果をリストで返すという処理のために、リスト内包表記というのがありますね。
 ある教科書に載ってた簡単な例でいうと、

cities = ["OSAKA", "TOKYO", "KYOTO"]


 というリストがあったとして、これを全部小文字に変えたいとします。
 そういう時、

cities_new = []
for city in cities:
	city = city.lower()
	cities_new.append(city)


 みたいなことをやらずに(やってもいいが)、

cities_new = [city.lower() for city in cities]


 と書けば終了というやつです。
 実際やってみると、

>>> cities_new = [city.lower() for city in cities]
>>> print(cities_new)
['osaka', 'tokyo', 'kyoto']


 ちゃんとできてますね。
 意味合いとしては、citiesというリストに入っている(in cities)個々の要素にcityという仮の名前をつけて、このcityに対して(for city)、順番にcity.lower()という処理を施して、返される結果が並んだリストを得る、ということですね。


 ちなみに、in citiesのあとにif文をつなげて、条件付きにすることもできます。

>>> cities_new = [city.lower() for city in cities if city.startswith("T")]
>>> print(cities_new)
['tokyo']

 
 
 ところでRでも、リスト*1やベクトル*2に入っている要素に対して、順番に何かの処理を施して、その結果をリストやベクトルで得たいときはあります。
 それでオライリーの『入門機械学習』という、機械学習の基本的な処理をRで演習する教科書を以前読んだのですが、その中ではapplyファミリーを使った↓のような書き方が多用されてました。
 上記の例でやるとすると、

cities <- c("OSAKA", "TOKYO", "KYOTO")
cities_new <- sapply(cities, function(p) {tolower(p)})


 というふうに書きます。ベクトルに対してつかうときはsapply()、リストに対してつかうときはlapply()になります。
 sapply()やlapply()は、1つめの引数で与えたベクトルやリストの要素1個1個について、2つめの引数で与えた関数を適用した結果をまとめて返します。で、この第2引数のところにpという仮の引数を持つ関数定義を埋め込んでやれば、既存の関数を使わなくても色々な処理を自由に書いて、リスト内包表記みたいな処理を一発で書くことができるわけですね。


 上記コードをRで動かしてみたら、

> print(cities_new)
  OSAKA   TOKYO   KYOTO 
"osaka" "tokyo" "kyoto" 


 なんか要素名称もついてきましたが・・・まあやりたいことはできています。
 
 

2.lm()やglm()を内包する関数定義をするとき

 lmでもglmでもいいんですが、関数定義の中に埋め込みたい場合があります。
 たとえば回帰モデルを作って当てはめを行い、その結果を報告しやすい形にまとめて返す関数を作っておきたいとかですね。具体的には、データフレームと、目的変数の列名称、説明変数の列名称を引数として与えると、重回帰分析を行って、summary(lm())で出てくる結果の表の一部を取り出してCSVで出力する・・・みたいな。
 ちなみに私は、たとえばp値に対応する記号を、lm()のデフォルトの"***", "**", "*", "."ではなく、"**", "*", "†"に書き換える、みたいなことを一括でやりたい場合とかがありました(†はascii文字ではないので文字コードの問題に注意が必要ですが)。


 以下では面倒なので、単純にsummary(lm())で出てくる表を表示するところまででやってみますが、そういう関数定義をやろうとした場合に、たとえば、

get_lmsummary1 <- function(df, res, v1, v2) {
	model <- lm(formula=res ~ v1 + v2, data=df)
	result <- summary(model)$coefficients
	return(result)
}


 みたいな書き方をしてもきちんと動きません。以下実際に試してみます。


 まず実験用のデータを乱数で適当につくります。

> set.seed(1)
> predictor_1 = rnorm(n=100, mean=5, sd=3)
> predictor_2 = rnorm(n=100, mean=10, sd=5)
> response =  3.0 + (1.5 * predictor_1) + (2.0 * predictor_2) + rnorm(n=100, mean=0, sd=10)
> 
> d <- data.frame(p_1 = predictor_1, 
+                 p_2 = predictor_2, 
+                 r    = response)
> 
> head(d)
       p_1       p_2         r
1 3.120639  6.898167 25.571309
2 5.550930 10.210579 48.636287
3 2.493114  5.445392 33.496339
4 9.785842 10.790144 35.949973
5 5.988523  6.727077  2.584583
6 2.538595 18.836436 69.457381


 一応、ためしにふつうに回帰分析してみます。

> model <- lm(formula=r ~ p_1 + p_2 , data=d)
> summary(model)

Call:
lm(formula = r ~ p_1 + p_2, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.4359  -4.3645   0.0202   6.3692  26.3941 

Coefficients:
            Estimate Std. Error t value          Pr(>|t|)    
(Intercept)   3.9710     3.1638   1.255          0.212440    
p_1           1.5704     0.3892   4.035          0.000109 ***
p_2           1.8931     0.2190   8.646 0.000000000000112 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.43 on 97 degrees of freedom
Multiple R-squared:  0.4839,	Adjusted R-squared:  0.4733 
F-statistic: 45.48 on 2 and 97 DF,  p-value: 0.00000000000001164


 さてこのデータフレームにさっきの関数を適用すると、
 

> get_lmsummary1(df=d, res=r, v1=p_1, v2=p_2)
Error in eval(expr, envir, enclos) : object 'r' not found


 rという名前のオブジェクトはない、とエラーがでます。まぁ、なるほど。
 引数を"r"のように文字列で与えるようにしてもダメでした。


 で、なんかいい解決方法あるのかもしれませんが、よくわからないので、私はとりあえず変数名を文字列で与えて、dfの列名称を検索して列番号を得ることにしました。ムリヤリ感がかなりあります。

> get_lmsummary2 <- function(df, res, v1, v2) {
+ 	res_n <- grep(res, colnames(df))
+ 	v1_n <- grep(v1, colnames(df))
+ 	v2_n <- grep(v2, colnames(df))
+ 	model <- lm(formula=df[,res_n] ~ df[,v1_n] + df[,v2_n])
+ 	result <- summary(model)$coefficients
+ 	rownames(result)[2:3] <- c(colnames(df)[v1_n], colnames(df)[v2_n])
+ 	return(result)
+ }
> 
> get_lmsummary2(df=d, res="r", v1="p_1", v2="p_2")
            Estimate Std. Error  t value           Pr(>|t|)
(Intercept) 3.971034  3.1637959 1.255149 0.2124402378606492
p_1         1.570367  0.3891757 4.035111 0.0001089365648451
p_2         1.893066  0.2189571 8.645832 0.0000000000001118


 これで思った通りの動きにはなりました。


 注意点として、grep()で当てているので、似たような列名称がある場合(1つの列名称をと同じ文字列を含む別の列が存在する場合)に、エラーが起きます。
 たとえば2個目の説明変数の列名称を"p_1_2"にしてみると、"p_1"というパターンは"p_1"にも"p_1_2"にもマッチするので、困ったことになります。

> d2 <- data.frame(p_1 = predictor_1, 
+                 p_1_2 = predictor_2, 
+                 r    = response)
> 
> get_lmsummary2(df=d2, res="r", v1="p_1", v2="p_1_2")
Error in model.frame.default(formula = df[, res_n] ~ df[, v1_n] + df[,  : 
  invalid type (list) for variable 'df[, v1_n]'


 幸い、grep()は正規表現が使えるので、"^"と"$"で挟んで完全一致にしてやると、

> get_lmsummary2(df=d2, res="r", v1="^p_1$", v2="p_1_2")
            Estimate Std. Error  t value           Pr(>|t|)
(Intercept) 3.971034  3.1637959 1.255149 0.2124402378606492
p_1         1.570367  0.3891757 4.035111 0.0001089365648451
p_1_2       1.893066  0.2189571 8.645832 0.0000000000001118


 というふうにちゃんと動きます。
 とりあえず、何かムダなことをやってる感があるので、もうちょっと綺麗な方法があれば知りたいです。
 

*1:Pythonのリストとは意味合いが違う

*2:Pythonには組み込みでベクトルという型はない