Rのlmで決定係数が切片”あり”・”なし”で異なる件について調べてみた

RとExcelのグラフで同じ回帰式でも決定係数 $R^ 2$ が異なるが、自分はてっきりExcel側のバグかと思っていたので、正直気にしていなかったが、決定係数 R2の違い: Excel, OpenOffice, LibreOffice および統計解析ソフト R を用いてによるとソフトウェアごとの $R^ 2$ の定義が異なっているからとのこと。

そして参考文献（Kvalseth T. O. (1985): Cautionary Note about R2, The American Statistician 39: 279-285.）いわく、決定係数 $R^ 2$ には8種類異なる定義があるよう（添字をつけて $R_1^ 2$ ~ $R_8^ 2$ ）。

Excelがどういう処理をしているのかの検証は今回別件で、Rのlmの処理で切片の有無で決定係数 $R^ 2$ が異なる件について調べてみたので今の処分った範囲でメモとして。

もしかしたらどこかにきれいにドキュメントがまとめてあったのかもしれないけど、見つけられず。

summary.lmのR document曰く

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 (1985) でいうところの

$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}$

とのこと（数式表記はKvalseth (1985) にならって）。

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

ただし、切片がないというのはintercept = 0ではないよう。

決定係数 R2の違い: Excel, OpenOffice, LibreOffice および統計解析ソフト R を用いてでも触れられているけど、切片が0になっているというのと切片を意図的に”なし”にするのかで計算結果が異なる様子（多分、リンク先の決定係数は入れ替わってる？）。

どんな動作をしているのかsummary.lmのソースコードを見に行くと

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

となっていて、切片が0かどうかを見ているのではなく、interceptという言葉があるかどうかを見てinterceptがあれば前述の $R_1^ 2$ で計算し、なければ（つまりlm(y~x1+x2+...-1)で原点を通るように指定していると） $R_7^ 2$ で計算するということになっている。

~~なんでlmが条件付けて決定係数 $R^ 2$ の計算を変えているのかはわからなかった。~~

ちょっと調べてもわからない感じだった。

多分こういう理由だろう、というのはあるけれども今回は調べてわかった範囲のみでメモとしておく。

~~今後の課題（いつか調べる？）。~~

（追記）CRANのFAQに理由とおぼしき記載があったが……

理由について言及されたCARNのFAQについて書いた。

indenkun.hatenablog.com

以下、検証用に作ったKvalseth (1985)の式を再現するようにしたポンコツ自作関数

$R_1^ 2$ を求める自作関数

lm.R21 <- function(form){
  z <- lm(form)
  RSS <- sum(z$residuals^2)
  MSS <- sum((z$fitted.values - mean(z$fitted.values))^2)
  TSS <- RSS + MSS
  1-(RSS/TSS)
}

使い方はlm.R21(y~x1+x2+...)とすると $R_1^ 2$ が出力される。

$R_7^ 2$ を求める自作関数

lm.R27 <- function(form){
  z <- lm(form)
  RSS <- sum(z$residuals^2)
  MSS <- sum(z$fitted.values^2)
  TSS <- RSS + MSS
  1-(RSS/TSS)
}

使い方はlm.R27(y~x1+x2+...)とすると $R_7^ 2$ が出力される。

備忘ログ

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

Rのlmで決定係数が切片”あり”・”なし”で異なる件について調べてみた