Rのlmの決定係数R^2が実数範囲で負をとりえない話

indenkun.hatenablog.com

の続き

Rのlmのsummary.lmの決定係数 $R^ 2$ についてながながと調べてきて、切片なしで計算式が異なること、その理由が決定係数 $R^ 2$ が負の値をとり得るからというところまで今までの流れ。

しかし、summary.lmのソースコードを読むと実数範囲ではRの計算方法だと負の値を取り得ないよねという話をしていく。

Rによるの計算をソースコードからみる
決定係数だと負の値をとり得る
決定係数だと負の値はとりえない
切片なしの回帰分析でRのsummary.lmにかかれている決定係数の式と実際にRで出力される決定係数が一致する近づく条件
ここまでわかったことまとめ

Rによる $R^ 2$ の計算をソースコードからみる

今まで調べたところによるとsummary.lmのドキュメントによると決定係数 $R^ 2$ は切片の有無で

r.squared
R², the ‘fraction of variance explained by the model’,

R² = 1 - Sum(R[i]^2) / Sum((y[i]- y*)²),

where y* is the mean of y[i] if there is an intercept and zero otherwise.

つまり、Rで決定係数 $R^ 2$ をsummary.lmで計算するとKvalseth T. O. (1985): Cautionary Note about R2, The American Statistician 39: 279-285.でいうところの

$R_1^ 2 = 1 - \dfrac{ \sum(y - \hat{y})^ 2 }{ \sum(y - \bar{y})^ 2}$

$R_7^ 2 = 1 - \dfrac{ \sum(y - \hat{y})^ 2 }{ \sum y^ 2}$

切片がある場合は、決定係数が $R_1^ 2$ で計算されるけど、切片なしの場合は $R_7^ 2$ で計算されるとのこと。

実際のRでの挙動としてはソースコードを見ると

r.squaredが

ans$r.squared <- mss/(mss + rss)

https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/lm.R lm.R v3.3.0

で計算されており、

$R^ 2 = \dfrac {mss}{mss+rss}$

mss: model sum of squares（回帰平方和）

rss: residual sum of squares（残差平方和）

ということ。

で、mssとrssは

    r <- z$residuals
    f <- z$fitted.values
    w <- z$weights
    if (is.null(w)) {
        mss <- if (attr(z$terms, "intercept"))
            sum((f - mean(f))^2) else sum(f^2)
        rss <- sum(r^2)
    } else {
        mss <- if (attr(z$terms, "intercept")) {
            m <- sum(w * f /sum(w))
            sum(w * (f - m)^2)
        } else sum(w * f^2)
        rss <- sum(w * r^2)
        r <- sqrt(w) * r
    }

https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/lm.R#L298 lm.R v3.3.0

で定義されているとのこと。

mssもrssも平方和なので実数範囲では負の値を取り得ないので、

$mss \geqq 0 \\ rss \geqq 0 \\ mss \leqq mss + rss \\ 1 \geqq \dfrac {mss}{mss+rss} \geqq 0$

なので、Rのlmの決定係数 $R^ 2$ は負の値を取り得ないとなる。

さらに言えば、Excelでたまに出る絶対値1を大きく上回る決定係数 $R^ 2$ は出力し得ない。

決定係数 $R_1^ 2$ だと負の値をとり得る

じゃあ、

$R_1^ 2 = 1 - \dfrac{ \sum(y - \hat{y})^ 2 }{ \sum(y - \bar{y})^ 2}$

だと決定係数 $R^ 2$ はどうなるのかというと、平方和は実数範囲では負の値を取りえないのは同じだけれども

${ \sum(y - \hat{y})^ 2 } > { \sum(y - \bar{y})^ 2}$

の条件さえ満たせば容易に負の値をとり得る。

この条件当てはまりが悪い（残差平方和が大きい）回帰分析を行うと、平均値の条件によっては容易に満たし得る。

それは今までの流れをくんで言うと、無理やり原点を通して回帰分析を行っている場合そうなり得るということ。

2変量の線形回帰分析についてRで $R_1^ 2 = 1 - \dfrac{ \sum(y - \hat{y})^ 2 }{ \sum(y - \bar{y})^ 2}$ をちゃんと再現すると

lm.R21.kavl <- function(x, y, intercept = c("existing", "zero")){
  z <- if(intercept == "existing")lm(y~x)
  else lm(y~x-1)
  R <- sum(z$residuals^2)
  1-(R/sum((y- mean(y))^2))
}

これで、lm(x,y,intercept = "zero")と入力するとlm(y~x-1)の原点を通る線形回帰分析の決定係数を $R_1^ 2$ で計算できる。

（”zero”って指定しないようにしても良かったけど、切片を意識したかったから）

ちなみにlm(x,y,intercept = "existing")と入力（intercept = "existing"はなくても計算できるが最後にちょっとしたエラーが出る）するとlm(y~x)の線形回帰分析の決定係数を $R_1^ 2$ で計算できる。

これを、決定係数 R2の違い: Excel, OpenOffice, LibreOffice および統計解析ソフト R を用いてで示された、負の値がでる例の

x <-  c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <-  c(1, 2, 2, 4, 3, 3, 2, 2, 3, 4)

のデータで実行するとintercept = "zero"とするとちゃんと

 lm.R21.kavl(x,y,"zero")
[1] -0.4257885

とexcelと同じ負の値がでるし、

x<- c(110, 120, 130, 140, 150, 160, 170, 180, 190, 200)
y<- c(180, 170, 180, 170, 160, 160, 150, 145, 140, 145)
lm.R21.kavl(x,y,"zero")
[1] -8.218253

とexcelと同じ絶対値1と大きく超える負の値がでる。

決定係数 $R_7^ 2$ だと負の値はとりえない

で、決定係数 $R_7^ 2$ だと実数範囲内では

$\sum(y - \hat{y})^ 2 \leqq \sum y^ 2$

になるので、負の値を取り得ないし、 $1 \geqq \dfrac{ \sum(y - \hat{y})^ 2 }{ \sum y^ 2} \geqq 0$ なので $1 \geqq 1 - \dfrac{ \sum(y - \hat{y})^ 2 }{ \sum y^ 2} \geqq 0$ となり絶対値1を超えることもないはず。

2変量の線形回帰分析についてRで $R_7^ 2 = 1 - \dfrac{ \sum(y - \hat{y})^ 2 }{ \sum y^ 2}$ をちゃんと再現すると

lm.R27.kavl <- function(x, y, intercept = c("existing", "zero")){
  z <- if(intercept == "existing")lm(y~x)
  else lm(y~x-1)
  R <- sum(z$residuals^2)
  1-(R/sum(y^2))
}

で上の関数と同じように使うと

x <- c(110, 120, 130, 140, 150, 160, 170, 180, 190, 200)
y <- c(180, 170, 180, 170, 160, 160, 150, 145, 140, 145)
lm.R27.kavl(x,y,"zero")
[1] 0.9303137

と $R_1^ 2$ では-8.218だったのに、 $R_7^ 2$ だと激アツな決定係数 $R^ 2$ が得られる。

$R_1^ 2$ が絶対値1を超える負の値を取るのに比較して、 $R_7^ 2$ だと絶対値1を超えない範囲で正の値を取るので一見妥当な決定係数 $R^ 2$ に見えるが、このx,y,をプロットして、原点を通過する回帰直線を描くと

plot(x,y)
abline(lm(x~y-1))

f:id:indenkun:20200304182202p:plain — fig.1

ということで全然、当てはまっていない。

R FAQやFAQ: Why are R2 and F so large for models without a constant?が触れているが、これは、そもそも決定係数の定義式の問題じゃなくx=0でy=0が想定されるものを無理やり原点を通る回帰式に当てはめようとしているせいだと思う。

なので、原点を通る回帰式でどの決定係数 $R^ 2$ の定義式が妥当かは（数学的な統計学上の問題はさておいて）、もともとのデータからあるいは経験則的に0点を通ることが妥当でない場合は、実務的にはそもそも妥当な回帰直線は得られていないんだから決定係数 $R^ 2$ がどうとかいう問題じゃないということになる。

切片なしの回帰分析でRのsummary.lmにかかれている決定係数 $R^ 2$ の式と実際にRで出力される決定係数 $R^ 2$ が一致する近づく条件

上記から推定されることとして、切片なしの回帰分析でRのsummary.lmにかかれている決定係数 $R^ 2$ の式と実際にRで出力される決定係数 $R^ 2$ が~~一致する~~近づく最低条件としては

（一致する条件はyの平均値が0であることが必要か？）

$\dfrac{ \sum(y - \hat{y})^ 2 } { \sum{y}^ 2} \leqq 1$

と考えられる。

この条件は残差平方和が $sum{y}^ 2$ よりも小さい（または同じ）という条件で、この条件を満たしているということは、どういう基準かにもよるけど原点を明らかに通らないデータで回帰直線を原点通してやっていいるよりは当てはまりが悪くがない可能性がある。

一方で

$\dfrac{ \sum(y - \hat{y})^ 2 } { \sum{y}^ 2} > 1$

この条件ではsummary.lmで定義されている決定係数[tex:R² ]と乖離してしまうが、残差平方和が $\sum{y}^ 2$ よりも大きいということで、どういう基準かにもよるがかなり当てはまりが悪い。つまり決定係数 $R^ 2$ という問題以前の回帰式が得られている可能性があるということになる。

（追記）同じデータでそれぞれ計算したときにどれくらいずれるかプロットしてみた。

indenkun.hatenablog.com

プロットしてみてわかったけど、 $R_1^ 2$ と $R_7^ 2$ はそもそも式が異なるので結構ばらつくパターンが多いと思った。

ここまでわかったことまとめ

決定係数 $R^ 2$ は複数の定義式がある 2.Rの決定係数 $R^ 2$ の定義式はsummary.lmのドキュメントに書いてある 3.Rの決定係数 $R^ 2$ は定義式とは（若干）異なる挙動で算出されている 4.Rの決定係数 $R^ 2$ は原点を通る回帰直線では当てはまりが著しく悪い場合summary.lmのドキュメントに記載されている式で算出される可能性がある

備忘ログ

チラシの裏的備忘録&メモ

Rのlmの決定係数R^2が実数範囲で負をとりえない話

Rによる $R^ 2$ の計算をソースコードからみる

決定係数 $R_1^ 2$ だと負の値をとり得る

決定係数 $R_7^ 2$ だと負の値はとりえない

切片なしの回帰分析でRのsummary.lmにかかれている決定係数 $R^ 2$ の式と実際にRで出力される決定係数 $R^ 2$ が一致する近づく条件

ここまでわかったことまとめ

Rによるの計算をソースコードからみる

決定係数だと負の値をとり得る

決定係数だと負の値はとりえない

切片なしの回帰分析でRのsummary.lmにかかれている決定係数の式と実際にRで出力される決定係数が一致する近づく条件

ここまでわかったことまとめ

Rによる $R^ 2$ の計算をソースコードからみる

決定係数 $R_1^ 2$ だと負の値をとり得る

決定係数 $R_7^ 2$ だと負の値はとりえない

切片なしの回帰分析でRのsummary.lmにかかれている決定係数 $R^ 2$ の式と実際にRで出力される決定係数 $R^ 2$ が一致する近づく条件