Rを使ってきて発見した小ネタを3つ紹介。もしかしたらもっといいやり方があるかもしれませんが、知っていたら教えて頂けると助かります。
利用するデータはいつもの以下のニコニコ動画から取得したデータです。
nico1
> nico <- read.table("http://am-yu.net/wp-content/uploads/2013/11/nico1.txt", sep="\t", header=T)
> head(nico)
category rank id view comment mylist
1 エンターテイメント 1 sm2057168 11278314 6249436 138031
2 エンターテイメント 2 sm1857897 3987887 829307 82675
3 エンターテイメント 3 sm1089673 3768025 389382 37145
4 エンターテイメント 4 sm4660488 3713982 225500 72711
5 エンターテイメント 5 sm9342638 3303266 92691 2049
6 エンターテイメント 6 sm12584996 3213439 224000 65855
一つの値しかない列を抽出
Rで重回帰分析しようとしたら、”対比は 2 つ以上の水準を持つ因子にしか適用できません”というエラーが表示されることがありました。原因は単純で、質的変数(カテゴリー変数)として利用しようとしていた変数が、1つの値しかもっていないというのが原因だったようです(1だけしかないとか0だけしかないなど、数値のみの場合はこのようなエラーはでませんでした)。
以下は、カテゴリーがVOCALOIDのデータのみを抜き出し、それを再生数を目的変数として、コメント数とマイリスト数とカテゴリを説明変数としようとして得られた結果です。
> nico.vocaloid <- nico[nico$category=='VOCALOID',] > nico.vocaloid.lm <- lm(view~comment+mylist+category, data=nico.vocaloid) Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 対比は 2 つ以上の水準を持つ因子にしか適用できません
この例の場合だと、質的変数はカテゴリーだけなので、これが原因ということは分かるのですが、質的変数が複数あった場合にいったいどれがエラーの原因になっているのか、判断できません。
で、いろいろ調べてみるとfactorという関数でカテゴリーを要素としたベクトルを作成することができると分かり、さらにlevels関数でグループ化された値のみ取得できるということがわかりました(参考：R-Source)。
それを参考に作ったのが以下。
> OnlyVariablePrint <- function(table){
+
+ len <- length(names(table))
+ for(VarCount in 1:len){
+ var <- length(levels(factor(table[,names(table)[VarCount]])))
+ if( var == 1 ){
+ print( names(table)[VarCount])
+ }
+ }
+ }
> OnlyVariablePrint(nico.vocaloid)
[1] "category"
変数が一種類しかない列があればその列名を表示するOnlyVariablePrintという関数を作成し、カテゴリーにVOCALOIDしか入ってないデータフレームを引数にいれると、"category"という出力が得られました。
やっていることは、"length(levels(factor(table[,names(table)[VarCount]])))"という箇所で、変数が何種類あるかどうか取得し、その後1つであればprint関数で名前を表示するという動きをしています。
モデル式から特定の変数を除外する
職場で扱っているデータは、かなりの数の変数が用意できるのですが、その中からできるかぎり精度の高い組み合わせというものを見つけなければいけません。とりあえず自分は、t値をみていって最も絶対値が小さいものを除外し、そのモデル式で再度重回帰分析をして評価し、先ほどより自由度調整済み決定係数が高ければ再度t値が最も小さいものを除外して再度重回帰分析を行い、先ほどより自由度調整済み決定係数が小さければ先ほど重回帰分析を行ったモデル式が最も精度が高いと判断するようなアルゴリズムを作成しました(調べてもその方法が正しいのかどうか怪しいところなのですが、今のところ、これで問題なさそうです)。
この時、モデル式から変数を除外するという作業を行わなければいけないわけですが、最初これをどうやればいいのか分かりませんでした。そこで、ふと、モデル式って文字列にしてみてもうまくいくのだろうかと思い、試してみるとできたのでなんとか解決することができました。
> ReduceFormula <- function(formula, CoefName){
+
+ formula2 <- formula
+ pasteCoef <- paste("[~]", CoefName, "[+]", sep="")
+ formula2 <- sub(pasteCoef, "~", formula2)
+ if(formula == formula2){
+ pasteCoef <- paste("[+]", CoefName, sep="")
+ formula2 <- sub(pasteCoef, "", formula2)
+ }
+ return(formula2)
+ }
> nico.formula <- "view~comment+mylist"
> ReduceFormula(nico.formula, "comment")
[1] "view~mylist"
> ReduceFormula(nico.formula, "mylist")
[1] "view~comment"
上記のReduceFormulaは、引数のモデル式から指定された変数を除外したモデル式を返す関数です。今回の例では、commentとmylistを説明変数としたモデル式を作成し、commentかmylistのみを説明変数から削除するといったことを行っています。
subという関数が正規表現でテキストを置換する関数です(最初に見つけた一つのみです)。前方一致の可能性があるため、除外しようとしていた変数が"~"の後につく変数だった場合、"+"の後につく変数かどうかの判断はしないようにしています(最初、if文をいれていないプログラムを書いてしまったことで、大変なことになりました)。"+"を"[+]"としているのは正規表現で判断しているためです。最初、"+"としていてうまくいかないので困りはてました。
本来なら、+の後につく変数があるかを判断したい場合、『pasteCoef <- paste("[+]", CoefName,"[+]", sep="")』と『pasteCoef <- paste("[+]", CoefName,"$", sep="")』の２パターン用意したほうがいいとは思うのですが・・・。正規表現が得意な人だともっと簡潔に書けるかもしれません。
特定の交互作用項のみの重回帰分析結果を取得する
前に交互作用モデルを含んだ重回帰分析について説明しました(交互作用モデルと、最初から変数として用意した重回帰分析の違いを調べてみた | while(isプログラマ))。
この時、最初から変数同士を掛けあわせた変数を用意しておくといい、なぜなら交互作用項のみの重回帰分析はできないようだから。なんて書いてしましましたが、すみません。嘘です。ちゃんとできました。
> nico.interaction.lm <- lm(view~comment:mylist, data=nico)
> summary(nico.interaction.lm)
Call:
lm(formula = view ~ comment:mylist, data = nico)
Residuals:
Min 1Q Median 3Q Max
-8300112 -368290 -197963 203775 7066614
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.335e+05 6.576e+03 81.13 <2e-16 ***
comment:mylist 1.190e-05 2.463e-07 48.31 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 622900 on 8998 degrees of freedom
Multiple R-squared: 0.206, Adjusted R-squared: 0.2059
F-statistic: 2334 on 1 and 8998 DF, p-value: < 2.2e-16
"*"ではなく、":"という記号を使うことでできるようです。"*"を使った時も結果は変数同士の間に":"があるからさっさと気づけ。という話です・・・。
ただ、これには少し問題があります。例えば『マイリスト数』と『コメント数とマイリスト数の交互作用項』を説明変数とした重回帰分析をしようとした場合、下記のようになります。
> nico.interaction.lm <- lm(view~mylist+comment:mylist, data=nico)
> summary(nico.interaction.lm)
Call:
lm(formula = view ~ mylist + comment:mylist, data = nico)
Residuals:
Min 1Q Median 3Q Max
-5397488 -205891 -82770 82872 4654344
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.379e+05 4.935e+03 48.20 <2e-16 ***
mylist 2.713e+01 2.376e-01 114.19 <2e-16 ***
mylist:comment 3.488e-06 1.737e-07 20.08 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 398000 on 8997 degrees of freedom
Multiple R-squared: 0.6758, Adjusted R-squared: 0.6757
F-statistic: 9378 on 2 and 8997 DF, p-value: < 2.2e-16
このように、モデル式のほうでは"comment:mylist"としているにも関わらず、結果のほうでは"mylist:comment"となってしまうのです。どうやら、単一の変数が説明変数に入っていれば、その変数が優先して先に書かれてしまうようです。
結果は同じなのだから別にいいじゃないか。と思われるかもしれませんが、そういうわけにはいきません。なぜなら、t値の絶対値が最も低い変数をモデル式から除外したい場合、モデル式のほうでは"comment:mylist"となっているのに、分析結果のほうでは"mylist:comment"となっているとt値がうまく取得できないためです。
> summary(nico.interaction.lm)$coefficients["mylist:comment","t value"] [1] 20.07621 > summary(nico.interaction.lm)$coefficients["comment:mylist","t value"] Error in summary(nico.interaction.lm)$coefficients["comment:mylist", "t value"] : subscript out of bounds
仕方がないので、片方でうまくいかない場合は、:の左右を交換して再度t値を取得するということをやっているのだけれども、もっとうまくやる方法はないものかどうか・・・。
やっぱり、それぞれの値を掛けあわせた変数を用意しておくのがいいのだろうか・・・。
コメント