Rの小技
Rでデータ分析するときに最近よくやる書き方があるのですが、よくやるといいながら1か月とか間が開くと忘れているので、メモしておきます。
２個あるのですが、１個目は将来ネットで検索して役に立ててくれる人がいるかもしれないので、１個目の小技をこのエントリのタイトルにして２個目をおまけ扱いにしました。
私はべつに、Rにもプログラミングにも詳しい人間ではないので、妥当なやり方かどうかは知りません。
というか、もっと良いやり方あればご教示いただけると・・・。（とくに２個目）
１．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"
なんか要素名称もついてきましたが・・・まあやりたいことはできています。
２．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()で当てているので、似たような列名称がある場合（１つの列名称をと同じ文字列を含む別の列が存在する場合）に、エラーが起きます。
たとえば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
というふうにちゃんと動きます。
とりあえず、何かムダなことをやってる感があるので、もうちょっと綺麗な方法があれば知りたいです。
