Rのsummary関数を使わずに、決定係数を求めてみる


前から言っているように仕事の業務で、重回帰分析を用いた解析を行っているのですが、ときどき予測値がありえない値になってしまうことがあります。具体的にいうと、マイナスの値にはならないのに、予測値がマイナスの数値をだしてしまうということです。
で、上からマイナスの場合は特定の値、例えば0にするようにということでした(実際にはもっと複雑なことを言われたのですが、ここでは単純に0とします)。まあ、それは想定内の指示だったのですが、ちょっと想定外だったのが、解析したときに予測値がマイナスの値になったものを0にして自由度調整済み決定係数を算出してほしいとのことでした。
普段、自由度調整済み決定係数はsummary関数で求めているわけですが、予測値を変更するとなると自分で求めなければいけません。
というわけで、今回は決定係数、自由度調整済み決定係数の求め方について勉強したことを書いていこうと思います。
使うデータは相変わらず、下記のニコニコ動画から取得したデータ。
nico

まずは、前準備。

データを取得して、解析しました。普段、コメント数とマイリスト数から再生数を予測するような解析をしていますが、今回はコメント数と再生数からマイリスト数を予測する解析をしていくことにします。
決定係数は0.6638783、自由度調整済み決定係数は0.6638036となりました。予測値をsummary関数でみてみると、どうやら値の低い予測値はマイナスの値になっているようです。
それでは、以下、自力で二つの決定係数を求めることに。

決定係数

決定係数 – Wikipedia
決定係数と自由度調整済み決定係数の求め方は、上記Wikipediaのページに書いてあります。これを参考に算出してみることに。
どうやら、1から、『残差二乗和』から『実測値と実測値の平均の二乗の和』を割った値を引いたものが決定係数らしい。
残差と実測値の平均を求めることができたら簡単に計算できそうなので下記のようにやってみました。

0.6638783とでました。この計算方法であっているようです。

自由度調整済み決定係数

自由度調整済み決定係数の求め方も先ほど記したWikipediaのページに書いてあります。
ただ、正直ぱっと見何をやっているのかよく分かりませんので、もう少し整理してみます。
ここで仮に説明変数の数をp、標本数の数をn、決定係数をR2とすると、自由度調整済み決定係数は下記のようになります。

というわけで、自由度調整済み決定係数も求めてみます。なお、説明変数の数はlmオブジェクトのrankというラベルから1を引いたものにあたり、標本数はnrow関数で取得できるので、以下のようにすると自由度調整済み決定係数も求められます。

というわけで、自由度調整済み決定係数も求めることができました。

では、試しに、予測値がマイナスになっているものは0にして決定係数を求めてみようと思います。
まず、予測値がマイナスになっているものは0にしたベクトルを作成。

つづいて、実測値と上記で作成した予測値の二つを含んだデータフレームを作成。

5番目の予測値がおかしい気がしますが、他はまああってそうです。

では、決定係数を求めてみることに。

決定係数は0.664791となり、少し増えました。

ここまでくれば後は同じなのですが、一応自由度調整済み決定係数も求めると下記のようになりました。

もちろん、自由度調整済み決定係数も増えています。
はたしてこれは、自由度調整済み決定係数といってしまっていいのかという疑問はありますが・・・。

実は、重回帰分析のやり方についてもlm関数を使わないでできないものかと勉強がてら調べてみたのですが、まだちょっとよく分かってません。想像以上に複雑でして・・・。
重回帰分析の仕組みをもっと分かるようになれば、何か新しい発見がある気もするのだけれども・・・。

Rの小ネタ3つ


Rを使ってきて発見した小ネタを3つ紹介。もしかしたらもっといいやり方があるかもしれませんが、知っていたら教えて頂けると助かります。

利用するデータはいつもの以下のニコニコ動画から取得したデータです。
nico1

一つの値しかない列を抽出

Rで重回帰分析しようとしたら、”対比は 2 つ以上の水準を持つ因子にしか適用できません”というエラーが表示されることがありました。原因は単純で、質的変数(カテゴリー変数)として利用しようとしていた変数が、1つの値しかもっていないというのが原因だったようです(1だけしかないとか0だけしかないなど、数値のみの場合はこのようなエラーはでませんでした)。
以下は、カテゴリーがVOCALOIDのデータのみを抜き出し、それを再生数を目的変数として、コメント数とマイリスト数とカテゴリを説明変数としようとして得られた結果です。

この例の場合だと、質的変数はカテゴリーだけなので、これが原因ということは分かるのですが、質的変数が複数あった場合にいったいどれがエラーの原因になっているのか、判断できません。
で、いろいろ調べてみるとfactorという関数でカテゴリーを要素としたベクトルを作成することができると分かり、さらにlevels関数でグループ化された値のみ取得できるということがわかりました(参考:R-Source)。
それを参考に作ったのが以下。

変数が一種類しかない列があればその列名を表示するOnlyVariablePrintという関数を作成し、カテゴリーにVOCALOIDしか入ってないデータフレームを引数にいれると、”category”という出力が得られました。
やっていることは、”length(levels(factor(table[,names(table)[VarCount]])))”という箇所で、変数が何種類あるかどうか取得し、その後1つであればprint関数で名前を表示するという動きをしています。

モデル式から特定の変数を除外する

職場で扱っているデータは、かなりの数の変数が用意できるのですが、その中からできるかぎり精度の高い組み合わせというものを見つけなければいけません。とりあえず自分は、t値をみていって最も絶対値が小さいものを除外し、そのモデル式で再度重回帰分析をして評価し、先ほどより自由度調整済み決定係数が高ければ再度t値が最も小さいものを除外して再度重回帰分析を行い、先ほどより自由度調整済み決定係数が小さければ先ほど重回帰分析を行ったモデル式が最も精度が高いと判断するようなアルゴリズムを作成しました(調べてもその方法が正しいのかどうか怪しいところなのですが、今のところ、これで問題なさそうです)。
この時、モデル式から変数を除外するという作業を行わなければいけないわけですが、最初これをどうやればいいのか分かりませんでした。そこで、ふと、モデル式って文字列にしてみてもうまくいくのだろうかと思い、試してみるとできたのでなんとか解決することができました。

上記のReduceFormulaは、引数のモデル式から指定された変数を除外したモデル式を返す関数です。今回の例では、commentとmylistを説明変数としたモデル式を作成し、commentかmylistのみを説明変数から削除するといったことを行っています。
subという関数が正規表現でテキストを置換する関数です(最初に見つけた一つのみです)。前方一致の可能性があるため、除外しようとしていた変数が”~”の後につく変数だった場合、”+”の後につく変数かどうかの判断はしないようにしています(最初、if文をいれていないプログラムを書いてしまったことで、大変なことになりました)。”+”を”[+]”としているのは正規表現で判断しているためです。最初、”+”としていてうまくいかないので困りはてました。
本来なら、+の後につく変数があるかを判断したい場合、『pasteCoef <- paste("[+]", CoefName,"[+]", sep="")』と『pasteCoef <- paste("[+]", CoefName,"$", sep="")』の2パターン用意したほうがいいとは思うのですが・・・。正規表現が得意な人だともっと簡潔に書けるかもしれません。

特定の交互作用項のみの重回帰分析結果を取得する

前に交互作用モデルを含んだ重回帰分析について説明しました(交互作用モデルと、最初から変数として用意した重回帰分析の違いを調べてみた | while(isプログラマ))。
この時、最初から変数同士を掛けあわせた変数を用意しておくといい、なぜなら交互作用項のみの重回帰分析はできないようだから。なんて書いてしましましたが、すみません。嘘です。ちゃんとできました。

“*”ではなく、”:”という記号を使うことでできるようです。”*”を使った時も結果は変数同士の間に”:”があるからさっさと気づけ。という話です・・・。
ただ、これには少し問題があります。例えば『マイリスト数』と『コメント数とマイリスト数の交互作用項』を説明変数とした重回帰分析をしようとした場合、下記のようになります。

このように、モデル式のほうでは”comment:mylist”としているにも関わらず、結果のほうでは”mylist:comment”となってしまうのです。どうやら、単一の変数が説明変数に入っていれば、その変数が優先して先に書かれてしまうようです。
結果は同じなのだから別にいいじゃないか。と思われるかもしれませんが、そういうわけにはいきません。なぜなら、t値の絶対値が最も低い変数をモデル式から除外したい場合、モデル式のほうでは”comment:mylist”となっているのに、分析結果のほうでは”mylist:comment”となっているとt値がうまく取得できないためです。

仕方がないので、片方でうまくいかない場合は、:の左右を交換して再度t値を取得するということをやっているのだけれども、もっとうまくやる方法はないものかどうか・・・。
やっぱり、それぞれの値を掛けあわせた変数を用意しておくのがいいのだろうか・・・。

Rでウェブスクレイピングしてみた


Rでウェブスクレイピングがやってみたくなったので調べてみた。
データサイエンティスト養成読本』でも、ウェブスクレイピングの章はRではなく、Pythonでの解説だったので、もしかしてRだと難しいのか? と思ったけど、XMLでのデータ取得とほぼ同じようにできた。

やり方は前にRでニコニコ動画APIを取得してみる | while(isプログラマ)で説明したやり方とほぼ同じ。
というわけで、今回ははてなブックマークのニコニコ動画の新着エントリーから、タイトルとURLとブクマ数を取得してデータフレーム化するプログラムを書いてみた。

できたテキストファイルが以下。
nicohatebu

XMLの時との違いは、xmlParseではなくて、htmlParseぐらいで、他は同じなよう。
簡単に何をやっているか説明すると、タイトルはh3タグの中のaタグの中のテキストなので、それを6行目でそのように取得。URLはh3タグの中のaタグの中のhref属性なので、それを7行目で取得。ブクマ数はliタグのdata-bookmark-count属性に書いてあるのでそれを8行目で取得しています。
あまり安全ではない書き方なような気がしますが、今のところこれで大丈夫そうなのでとりあえず。
本当は、個別個別のliをまず取得して、その中からh3タグの中のaタグを取得。というようにしたかったのだけれども、やり方が分かりませんでした。下記のようにやるとできると思ったのですが、うまくいかなくて・・・。

[[1]]の”【疑似m@s】割と熱湯のお風呂 ‐ ニコニコ動画:GINZA”しか取得しないと想定してやってみたのですが・・・。ここでは、『li[[1]]』が外部ポインターオブジェクトというものらしいのですが、まだよく分からないでいます・・・。

交互作用モデルと、最初から変数として用意した重回帰分析の違いを調べてみた


まだ試してませんが、仕事でどうも変数同士が交互作用しているような気がする変数があります。例えば、中古の携帯電話の価格の推定をするとした場合、説明変数に発売年とカメラの画素数を加えるとすると、去年発売されたものの画素数が200万ぐらいなら価値がかなり下がるが、10年ぐらい前だとむしろ価値があがるといった具合です。でも、これを普通の回帰分析でやってしまうと、画素数の係数は一つの値しかでないので、発売年が去年でも10年前でも同じ値だけプラスすることになります。
で、これを解決してくれそうなのが交互作用モデルを用いた重回帰分析というわけです。
いや、実は言うと交互作用モデルについてはよく分かってない部分があるので、本当にこの解釈があっているかどうか自身を持ってはいえないのですが、そんなに間違ってないとは思います。多分。

ちなみに、交互作用モデルというのを簡単にいうと、上の携帯電話の価格の例でいうと、普通の重回帰分析では下記のような式を考えます。
『予測価格=a*発売年+b*カメラの画素数』
普通の重回帰分析ではこのaとbを求めるわけですが、交互作用モデルを用いた重回帰分析では下記のような式を考えます。
『予測価格=a*発売年+b*カメラの画素数+c*発売年*カメラの画素数』
この、aとbとcを求めるわけです。

ただ、ここで思ったわけです。最初から、『発売年*カメラの画素数』というような変数があった場合は、結果が異なるのかどうなのかと。例えば、『発売年*カメラの画素数』を含んだ『年画素』という変数が会った場合、
『予測価格=a’*発売年+b’*カメラの画素数+c’*年画素』
という式になるaとa’、bとb’、cとc’は違うのかどうなのかと。カテゴリー変数得られた係数の数値をそれぞれのカテゴリーに割り当てて重回帰分析した結果ですら同じになったので、多分同じになるんじゃないかと思いながらもやってみました。

利用したファイルは下記。
nico4
先ほど、ニコニコ動画の24時間の総合ランキングから取得したデータです。普段使っているニコ動のそれぞれのカテゴリーの合計データを使おうと思ったのですが、オーバーフローしてしまったので新しくとってきました。
今回試したのは、コメント数とマイリスト数を用いて再生数を推定するというもの。
結果は以下。

最初の『view~comment+mylist+comment*mylist』は『view~comment*mylist』でもOKです。
というわけで、結論。
全く同じなようです。自由度調整済み決定係数ぐらい異なるかもしれないと思っていましたが、それすら同じでした。

最初から変数同士の掛け算を入れた変数を用意できるなら、そうやるほうがいいような気がします。これだと、Rで、コメント数とコメント数とマイリスト数を掛けた値のみの重回帰分析というのも簡単にできるので(自分がやり方を知らないだけだとは思いますが)。

Rから回帰分析のいろいろな情報を取得する方法


Rを使って重回帰分析をして、summary関数で重回帰分析に関する情報が見れるとわかったものの、『係数だけ』『自由度調整済み決定係数だけ』がほしい。そう思って取得方法をいろいろ調べてみました。
例によって、使うのは下記のニコ動からのデータ。
nico

上記は、コメント数とマイリスト数とカテゴリーを説明変数にして、再生数を目的変数とした重回帰分析の結果。

係数

lmオブジェクトの後に『$coefficients』、もしくは『[“coefficients”]』、coef(coefficients)関数を用いても取得できる。

一つ一つを取得したい場合は、後ろに『[“”]』でくくって、中に変数名をいれる。

残差

lmオブジェクトに$residualをつけるか、residuals(resid)関数を利用

lmオブジェクトを引数にしたsummary関数で表示されるような四分位数を取得したい場合は、上記の記述をsummary関数の引数とする。

予測値

lmオブジェクトに$fitted.valuesをつけるか、predict関数の引数にlmオブジェクトをいれる

切片を含む、説明変数の数

lmオブジェクトに『$rank』

自由度(標本数から説明変数(切片含む)の数を引いたもの)

$df.residualをつける。

カテゴリー変数の値

$xlevelsをつける。

これを、”category”という文字列(もしくは、その文字列をいれた変数)から取得したいのだけどやりかたが分からない。nico.lm[“category”][1]では取得できないようだし。data.frame関数通したらできるようではあるのだけど。

まわりくどい。自分が知らないだけで、何かいいやり方があると思うのだけど。

モデル式

formula関数を使うか、$callをつける(結果は異なる)

AIC(赤池情報量規準)

AIC関数の引数にいれる

解析に使用したデータのみのモデル

$modelをつける。

決定係数

lmオブジェクトを引数にしたsummary関数の後ろに$r.squared

自由度調整済み決定係数

lmオブジェクトを引数にしたsummary関数の後ろに$adj.r.squared

標準誤差

lmオブジェクトを引数にしたsummary関数の後ろに$coefficients[,”Std. Error”]

t値

lmオブジェクトを引数にしたsummary関数の後ろに$coefficients[,”t value”]

p値

lmオブジェクトを引数にしたsummary関数の後ろに$coefficients[,”t value”]

他にもいろいろあるようですが、自分が必要になりそうなもののみ列挙しました。なお、ここに書いてあるのはネットや書籍で調べたわけではなく、names関数を使って調べて行った結果、分かった方法です。

names関数、すごい便利です。