前から言っているように仕事の業務で、重回帰分析を用いた解析を行っているのですが、ときどき予測値がありえない値になってしまうことがあります。具体的にいうと、マイナスの値にはならないのに、予測値がマイナスの数値をだしてしまうということです。
で、上からマイナスの場合は特定の値、例えば0にするようにということでした(実際にはもっと複雑なことを言われたのですが、ここでは単純に0とします)。まあ、それは想定内の指示だったのですが、ちょっと想定外だったのが、解析したときに予測値がマイナスの値になったものを0にして自由度調整済み決定係数を算出してほしいとのことでした。
普段、自由度調整済み決定係数はsummary関数で求めているわけですが、予測値を変更するとなると自分で求めなければいけません。
というわけで、今回は決定係数、自由度調整済み決定係数の求め方について勉強したことを書いていこうと思います。
使うデータは相変わらず、下記のニコニコ動画から取得したデータ。
nico
まずは、前準備。
> nico <- read.table("http://am-yu.net/wp-content/uploads/2013/11/nico1.txt", sep="\t", header=T)
> nico.lm <- lm(mylist~view+comment, data=nico)
> summary(nico.lm)
Call:
lm(formula = mylist ~ view + comment, data = nico)
Residuals:
Min 1Q Median 3Q Max
-105235 -3818 -266 1998 186886
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.204e+03 1.519e+02 -7.927 2.52e-15 ***
view 2.297e-02 1.741e-04 131.970 < 2e-16 ***
comment -2.344e-03 2.819e-04 -8.316 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 11300 on 8997 degrees of freedom
Multiple R-squared: 0.6639, Adjusted R-squared: 0.6638
F-statistic: 8885 on 2 and 8997 DF, p-value: < 2.2e-16
> summary(nico.lm)$r.squared
[1] 0.6638783
> summary(nico.lm)$adj.r.squared
[1] 0.6638036
> summary(predict(nico.lm))
Min. 1st Qu. Median Mean 3rd Qu. Max.
-27260 2564 6489 11350 15820 382100
データを取得して、解析しました。普段、コメント数とマイリスト数から再生数を予測するような解析をしていますが、今回はコメント数と再生数からマイリスト数を予測する解析をしていくことにします。
決定係数は0.6638783、自由度調整済み決定係数は0.6638036となりました。予測値をsummary関数でみてみると、どうやら値の低い予測値はマイナスの値になっているようです。
それでは、以下、自力で二つの決定係数を求めることに。
決定係数
決定係数 – Wikipedia
決定係数と自由度調整済み決定係数の求め方は、上記Wikipediaのページに書いてあります。これを参考に算出してみることに。
どうやら、1から、『残差二乗和』から『実測値と実測値の平均の二乗の和』を割った値を引いたものが決定係数らしい。
残差と実測値の平均を求めることができたら簡単に計算できそうなので下記のようにやってみました。
> nico.lm.rsquared <- 1-sum(residuals(nico.lm)^2)/sum((nico$mylist-mean(nico$mylist))^2) > nico.lm.rsquared [1] 0.6638783
0.6638783とでました。この計算方法であっているようです。
自由度調整済み決定係数
自由度調整済み決定係数の求め方も先ほど記したWikipediaのページに書いてあります。
ただ、正直ぱっと見何をやっているのかよく分かりませんので、もう少し整理してみます。
ここで仮に説明変数の数をp、標本数の数をn、決定係数をR2とすると、自由度調整済み決定係数は下記のようになります。
((n-1)*R2-p)/(n-1-p)
というわけで、自由度調整済み決定係数も求めてみます。なお、説明変数の数はlmオブジェクトのrankというラベルから1を引いたものにあたり、標本数はnrow関数で取得できるので、以下のようにすると自由度調整済み決定係数も求められます。
> nico.lm.adj.rsquared <- ((nrow(nico)-1)*nico.lm.rsquared-(nico.lm$rank-1))/((nrow(nico)-1)-(nico.lm$rank-1)) > nico.lm.adj.rsquared [1] 0.6638036
というわけで、自由度調整済み決定係数も求めることができました。
では、試しに、予測値がマイナスになっているものは0にして決定係数を求めてみようと思います。
まず、予測値がマイナスになっているものは0にしたベクトルを作成。
> nico.mylist.predict <- ifelse(predict(nico.lm)<0, 0, predict(nico.lm))
> summary(nico.mylist.predict)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0 2564 6489 11430 15820 382100
つづいて、実測値と上記で作成した予測値の二つを含んだデータフレームを作成。
> nico.dataframe <- data.frame(mylist=nico$mylist, predict=nico.mylist.predict) > head(nico.dataframe) mylist predict 1 138031 243265.56 2 82675 88473.21 3 37145 84453.03 4 72711 83595.51 5 2049 74470.67 6 65855 72099.15
5番目の予測値がおかしい気がしますが、他はまああってそうです。
では、決定係数を求めてみることに。
> nico.dataframe.rsquared <- 1-sum((nico.dataframe$mylist-nico.dataframe$predict)^2)/sum((nico.dataframe$mylist-mean(nico.dataframe$mylist))^2) > nico.dataframe.rsquared [1] 0.664791
決定係数は0.664791となり、少し増えました。
ここまでくれば後は同じなのですが、一応自由度調整済み決定係数も求めると下記のようになりました。
> nico.dataframe.adj.rsquared <- ((nrow(nico.dataframe)-1)*nico.dataframe.rsquared-(nico.lm$rank-1))/((nrow(nico.dataframe)-1)-(nico.lm$rank-1)) > nico.dataframe.adj.rsquared [1] 0.6647165
もちろん、自由度調整済み決定係数も増えています。
はたしてこれは、自由度調整済み決定係数といってしまっていいのかという疑問はありますが・・・。
実は、重回帰分析のやり方についてもlm関数を使わないでできないものかと勉強がてら調べてみたのですが、まだちょっとよく分かってません。想像以上に複雑でして・・・。
重回帰分析の仕組みをもっと分かるようになれば、何か新しい発見がある気もするのだけれども・・・。
コメント