2012年2月5日日曜日

地獄のミサワ的Rユーザーの回帰分析

このエントリーをはてなブックマークに追加
Clip to Evernote
Pocket

Rが使えるフリをするための14の知識」に対して、古参R使いから「Rユーザーの「フリ」なんてしなくていいんですよ」とコメントを頂いた。確かに、意識が高いRユーザーのフリをする人間は大変にウザイので、いない方がいいかも知れない。

どのぐらいウザイかは、身近にいるウザイ人物がRユーザーになったらどうなるかを想像すればすぐ分かる。例えば、そこそこ可愛い女子大生が、太めでマニアックでラーメンにこだわるR使いに質問すると、典型的には以下のような会話が展開されるであろう。

男「あれ、アイナちゃん、こんな所でどうしたの?もしかして、俺に告白したいの?」

女「・・・(キモヲタデブと付き合うぐらいなら、風俗で働くわボケ)。ミサワ君、こんにちは。」

男「何度も言っているけど、アイナちゃんの気持ちは嬉しいけど、全国のオレのファンが君を殺す事になるから分かってくれよ。」

女「・・・(どこにファンがいるんだよ、冬なのに変な汗が出てくる)。学校の課題のデータをRで分析しないといけないんだけど、ミサワ君、R得意だって言っていたよね?」

男「まぁね。StataでもSASも他の人の10倍ぐらいは得意だけど、Rは10年も使っているから100倍ぐらい得意かな。おいおい、日本のRの第一人者なんて言わないでくれよ、恥ずかしいから」

女「・・・(金積まれても言わねーーーーーーーーよ、ボケナス)。そ、その自信に満ち溢れすぎた知識で、重回帰を教えてくれない?」

男「教わろうとするな自分で学べが信条だけど、アイナちゃんの頼みとあっては取り下げるしかないな。OLSでいいのかな?」

女「・・・(教えたがりのアンタが他の女に数学を教えていたのは知っている)。OLS?わかんない。」

男「一番基本となる分析方法。オレみたいなエキスパートだと、まず使わないけどね。」

女「・・・(いちいち自分をエキスパート言うな!!!)。たぶん、それ。」

男「データの取り込みはできた?」

女「うん。read.table()でできた。」

男「さすがオレの女。」

女「・・・(いつから私はてめぇの女になったんだ?)。重回帰の仕方が分からない。」

男「OLSに限らず、まずはモデル式を立てる。アイナちゃんはモデルやったりするんだっけ?」

女「・・・(Rの話をしろよ)。どうやって立てるの?」

男「データフレームdframeのy、x1、x2がデータみたいだから、まずはdframeを省略するためにattach(dframe)としよう。次にfrml <- y ~ x1 + x2と書く。」

女「うんうん(顔、よせてくんな!)。」

男「次に説明変数を行列にしよう。X <- model.matrix(frml)。」

女「うんうん(あれ?授業で行列なんか使ったっけ?)。」

男「Xをprintすれば分かるけど、切片項も入っている。Intersectionね。」

女「うんうん(・・・難しいなぁ)」

男「行列で推定式を表すとY = Xβ + ε。これをX'Y = X'Xβ + X'ε、(X'X)-1X'Y = β + (X'X)-1X'εと整理できる。(X'X)-1X'εは平均0の誤差項だから、とりあえずは気にしない。」

女「う・・・」

男「Rで行列演算を書くと、b <- solve(t(X)%*%X)%*%t(X)%*%yとなり、bが係数になる。オレは線形代数は自家薬籠中の物だけど、アイナちゃんには難しいかな?」

女「・・・(じかやくろうちゅうって何?)。課題を出せばいいだけだから、頑張る。課題には係数の他に、標準偏差とt値とP値も出すって書いてあるのだけど・・・。」

男「係数以外も出すのか。しょーがないな、アイナちゃんのためだから頑張るよ。自由度dfはサンプル数から説明変数C、x1、x2の数3を引いた297ね。」

女「うん」

男「分散s2はsum((y - X %*% b)^2)/dfになるから、分散共分散行列vcovはs2*solve(t(X)%*%X)とかける」

女「うんうん(もう、計算できればいい)。」

男「標準偏差は分散共分散行列の対角成分のルートとったものだからsqrt(diag(vcov))ね。x1の係数のt値などを出すには、マスク行列tmをmatrix(c(0, 1, 0), 1, 3)としてまず作る。t値tvは(tm %*% b) / sqrt(tm %*% vcov %*% t(tm))で出せて、P値は両側検定としてif(0>tv) tp <- 2*pt(tv, df) else tp <- 2*(1 - pt(tv, df))。」

女「う・・・良く分からないけど、入力。」

男「P値は0に限りなく近い、つまりアイナちゃんとオレみたいな関係だから、sprintf("%f [%0.5f%%]", tv, tp/100)としないと0としか表示されないね」

女「・・・(近いのは今だけで、それも一生分としても十分だから)」

男「切片項とx2に関しても出してみてよ。」

女「お、出た!これ、計算あっているの?」

男「もう半年ぐらい計算を間違えたことはないよ。オレってちょっと緻密過ぎ?」

女「・・・(顔よせてくんな!厚かましすぎ!)。よし!重相関係数も出したい。」

男「残差平方和ssrをsum((y - X %*% b)^2)、yの分散var_yをsum((y - mean(y))^2)とすると、定義より決定係数は1 - ssr/var_y、自由度調整済み決定係数は1 - (ssr/df)/(var_y/(300-1))、相関係数はこれらのルート。超かんた~ん。」

先輩「二人とも何をやっているの?ちょっといい?」

女「サエコ先輩、こんにちは。」

先輩「さっきから見ていたのだけど、summary(lm(y ~ x1 + x2, data=dframe))では駄目なの?」

女「え・・・?それだけ?」

先輩「答えが出ればいいんでしょう?」

男「そういえば、lm()があったね。オレはいつも最新の計量分析手法を使っているから、ビルトインの分析関数の事なんてすっかり忘れていたよ(アハアハハ・・・)」

女・先輩「死ねよ、オマエ。」

もう意識の高いRユーザーのフリは、公害に近い外部不経済があるのがお分かりであろう。どちらかと言うとRに使われているユーザーは、簡単に統計解析を行いたいユーザーの阻害になってしまう(ミサワ君とサエコ先輩のソースコード)。そう、“Rユーザーの「フリ」なんてしなくていいんですよ”。

現実的にはデータの整理や解析結果の解釈に時間を使うべきで、統計手法のコーディングは気にしない方が望ましい。そういう意味ではRは実践の中で自然に覚えていくのが理想であり、無理して意識の高いユーザーを振舞うのは良い事では無いのは間違いない。

0 コメント:

コメントを投稿