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関数、すごい便利です。

自由度調整済み決定係数が小さくなるから、カテゴリー変数はダミー変数化しないほうがいい?


最近、Rについてのエントリーが多いことからも察している人も多いと思いますが、最近は仕事でRばっかり使ってます。ある製品群の値段予測システムという社長の案があって、自分はRを用いてその予測価格の算出方法をいろいろ考えているところです。そのデータ自体はデータベースに大量にあるので、それを有効活用しようというわけです。
前々から社長はその案を画策していたそうで、社長から渡された多変量解析の参考書籍には2011年に購入したであろうと思われる栞がはさんでありました。確かあちこちにメモ書きもあったと思うので、社長もよく勉強してるんだと思います。

ところで、解析の対象の項目には、順序があるわけではない、あるカテゴリー変数がありました。隠すことではないのと思うで言ってしまうと、製品の『色』です。ただ、だいたいの人気色というのは予想がつくので、自身でそれぞれの色に数値を割り当てていって精度をあげてほしいということでした。

試しに、下記のニコニコ動画から取得したデータを用いてやってみようと思います。
nico
この場合、ニコ動のカテゴリーがまんまカテゴリー変数ということなります(ややこしいですが)。まずは試しに、再生数とカテゴリーをboxplotで出力してみる。

すると、下記のような図となりました(y軸の制限をなくすと縦に長くなりすぎてわかりづらいので、箱が分かる範囲にしています)。
カテゴリーと再生数のボックスプロット
項目名が書いてないところが多いので分かりにくいですが、左から『R-18』『VOCALOID』『アイドルマスター』『アニメ』『エンターテイメント』『ゲーム』『スポーツ』『その他』『ニコニコインディーズ』『ニコニコ技術部』『ニコニコ手芸部』『ニコニコ動画講座』『ラジオ』『演奏してみた』『音楽』『科学』『歌ってみた』『作ってみた』『自然』『車載動画』『政治』『東方』『動物』『日記』『描いてみた』『踊ってみた』『旅行』『料理』『例のアレ』『歴史』 となっています(項目名を縦にして表示したいのですが、調べても分かりませんでした。改行を入れて対処するという方法もあるんですが、図とかぶってしまうので・・・)。
では、これを参考に数値化しています。とりあえず中心点を基準に上のほうにあるものを大きい数値、下の方にあるものを小さい数値であらわしてみます。

とりあえず下記のような感じにしてみました。

これで重回帰分析できます。試してみましょう。

決定係数(Multiple R-squared)が0.7105、自由度調整済み決定係数(Adjusted R-squared)が0.7104となっています。悪くなさそうです。もう少し微調整したらもっとよくなるかもしれません。

ところで、その時(一ヶ月ほど前ですが)は自分も統計に関しては全くもって無知で「確かに、数値化しないと計算できないしね。なるほどなるほど」なんて思ってやっていたのですが、いろいろ調べてみるとダミー変数化するという方法があるということを知りました。例えば、黒という変数を用意して、黒い製品には1、他は0。また、白という変数を用意して、白色なら1、ほかは0。という具合に、それぞれのカテゴリーごとに変数を作るというものです。これならわざわざ自分で順序を考えなくてもいいですし、楽ができそうです。

ということを社長に話してみると、こう言われました。
「変数が増えると補正R2が小さくなるからそれはやらないほうがいい」(補正R2とは、自由度調整済み決定係数のこと。Excelではこう記される。初めてExcelで調査してたのでこの言葉が社内では一般的)
なんとも腑に落ちない感じがしたのですが、自分も統計学に関しては初心者だったのでなんとも言えませんでした。
確かに、自由度調整済み決定係数は、決定係数や標本数の数が変わらない状況な場合、変数だけ増えると下がる傾向にあります。
ちなみに、変数の数をp、標本数の数をn、決定係数をR2とすると、自由度調整済み決定係数は下記のようになります。

参考:自由度調整済決定係数
上記の重回帰分析の結果を例にもとめてみると、((9000-1)*0.7105-3)/(9000-1-3)となり、確かに結果は0.7104となりました。

ではここで、標本数が20、R2が0.7、変数の数が9と10の場合を考えてみようと思います。
まず、変数の数が9の場合は((20-1)*0.7-10)/(20-1-10)となり、これは0.43になります。10の場合は0.3666666…となりました。確かに変数の数が増えると自由度調整済み決定件数は減るようです(確かあの時社長は、「分母が大きくなるから小さくなる」みたいなことを言っていたような気がするのですが・・・。記憶違いかもしれない。さすがに自分より一年以上も多く統計の勉強をしていた社長がそんな間違いをするとは思えないし)。

なお、上記のニコ動の例をダミー変数化すると下記のような結果に(指数表記をなくしてます)。

自由度調整済み決定係数が、増えてるという・・・。

さて、例がものすっごい悪かったのですが、実は業務で使っていた解析では数値を調整することにより、自由度調整済み決定係数がダミー変数化する場合より大きくなりました。
ところで、ダミー変数を使うということは変数の値は1か0かということになります。ということは、ダミー変数を使って出てきた係数を使って数値化し、1という数値のほうを係数と考えてもいいんじゃないか。と思ったわけです。
例えば、以下の様な感じ。

思った通り、係数は1となり、自由度調整済みも大きくなりました。もし、自由度調整済み係数をもとめられているけれども、ダミー変数化すると自由度調整済み係数が下がってしまう場合は、上記のようにして求めた自由度調整済み係数を報告すればいいかもしれませんね(自分はそうしてます)。

ということはですよ。ということは、今、上記の解析では『コメントの係数*コメント数+マイリストの係数*マイリスト数+ダミー変数の係数+切片』という式になっているわけですが、そもそも『コメントの係数*コメント数+マイリストの係数*マイリスト数+ダミー変数の係数』の部分を一つの変数と考えたらいいんじゃね? と。そしたら、変数の数は1つだけになり、自由度調整済み決定係数の値も大きくなるはず。
やってみましょう。

はい。確かに自由度調整済み決定係数が大きくなり、表示の上では決定係数と同じになりました(計算してみたら分かりますが、全く同じではありません。決定係数よりは小さくなります)。

自由度調整済み決定係数は違いますが、上記3つ(決定係数が0.7591となってる解析)は、どれも同じ予測値を算出します。だからといって参考にならないかというとそういうわけではないのですが、あまりアテにするのもよくないような気がします。少なくとも、サンプル数がかなり多ければ(別にどれぐらいの数以上かを想定してるわけではないです)、決定係数だけを見てもいいんじゃないかと。まあ、変数を一つ増やすか増やさないかの参考なら使えるかもしれませんが。

最後にもう一度言っておきますが、自分は統計学を勉強して1ヶ月程度しかたってない初心者の中の初心者です。もしかしたら盛大な勘違いをしているかもしれません。もし、おかしなことを書いてるなら指摘してください。