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

「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は実践の中で自然に覚えていくのが理想であり、無理して意識の高いユーザーを振舞うのは良い事では無いのは間違いない。