Rのlmの決定係数R^2について定義式の違いでの差がどの程度なのかプロットしてみた

indenkun.hatenablog.com

の続きの蛇足。

これまでの簡単なまとめ
Rのlmで求める式でのR21とR27の差について
Kvalsethの式で計算したR2の差
Kvalsethとlmの式でのRの差
結局何を言いたいのか

これまでの簡単なまとめ

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

Rのlmの決定係数R²（summary.lmで出力される結果）（tex表記で描くのがめんどくさいのでR²）は切片の有無で定義される式が異なります。

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.

（引用：summary.lm R Documentation）

これを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}$

となります。

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

以下tex表記がめんどくさいのでそれぞれR21、R27とします。

理由はCRANのFAQによると

7.41 Why does summary() report strange results for the R² estimate when I fit a linear model with no intercept?

As described in ?summary.lm, when the intercept is zero (e.g., from y ~ x - 1 or y ~ x + 0), summary.lm() uses the formula R² = 1 - Sum(R[i]^2) / Sum((y[i])²) which is different from the usual R² = 1 - Sum(R[i]^2) / Sum( (y[i] - mean(y) )²). There are several reasons for this:

Otherwise the R² could be negative (because the model with zero intercept can fit worse than the constant-mean model it is implicitly compared to).

If you set the slope to zero in the model with a line through the origin you get fitted values y*=0

The model with constant, non-zero mean is not nested in the model with a line through the origin. All these come down to saying that if you know a priori that E[Y]=0 when x=0 then the ‘null’ model that you should compare to the fitted line, the model where x doesn’t explain any of the variance, is the model where E[Y]=0 everywhere. (If you don’t know a priori that E[Y]=0 when x=0, then you probably shouldn’t be fitting a line through the origin.)

（引用：CRAN FAQ）

ということで、理由は決定係数R²が負の値を取らないようにとのことです。

しかし、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.0lm

ということでRのlmは定義式を直接計算しているわけではないので、Rのlmではどうしたっても負の値を取りえません。

この点については indenkun.hatenablog.com に書いた。

ここまでがこれまでの流れの簡単なまとめ。

そこで、R²について同じ2変数のxとyについて、Kvalsethの式にならったR21、R27とRのlmで計算しているR21、R27でどの程度差がでるのかシミュレーション（？）してみました。

KvalsethのR21とR27はそれぞれ、kavl_R21、kavl_27（名前の綴を間違っているのに今気づいたが、コードも全部それで書いたのでこのままいく）、Rのlmで計算されるR21、R27はそれぞれR_R21、R_R27とします。

（数学が得意だったらその差がどの条件でどの程度の範囲内に散布するのか計算できるかもしれないけど、そんな能力ないのでrnorm()で作った2変数についてループ処理で決定係数R²を求めていくスタイル。）

kvalsethの式を再現する自作関数は

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

で、x、yはデータで、R2=でR21、R27のどちらの式を抽出するか、intercept=で切片の有無（原点を通る回帰直線を強制的に書かせるか）を選択します。

Rのlmの式を再現する自作関数は

lm.R_R2 <- function(form, R2 = c("R21","R27")){
  z <- lm(form)
  RSS <- sum(z$residuals^2)
  if(R2 == "R21"){
    MSS <- sum((z$fitted.values - mean(z$fitted.values))^2)
  }
  else{
    MSS <- sum(z$fitted.values^2)
  }
  TSS <- RSS + MSS
  MSS/TSS
}

でfromでy~xとすると普通に回帰直線を、y~x-1などとすると原点を通る回帰直線について、R2=でR21、R27のどちらを求めるかを選択します。

Rのlmで求める式でのR21とR27の差について

計算用に使用した自作関数が、

lm.R2.RR <- function(n = 10, 
                     times = 100, 
                     intercept = c("existing", "zero")){
  R_R21 <- NULL
  R_R27 <- NULL
  if(intercept == "existing"){
    lmform <- y~x
  }
  else{
    lmform <- y~x-1
  }
  for(i in 1:times){
    x <- rnorm(n)
    y <- rnorm(n)
    R_R21[i] <- lm.R_R2(form = lmform, R2 = "R21")
    R_R27[i] <- lm.R_R2(form = lmform, R2 = "R27")
  }
  data.frame(R_R21, R_R27)
}

でRのlmで求められるR21とR27の差について2変量について10個ずつの乱数を発生させて、100回くらい（設定変えられる）で切片の有無で比較します。

今後plotに結果を直接渡すのにパイプ処理するので（単純に面倒なだけ）

library(magrittr)

でmagrittrのライブラリを読み込んでおいて、%>%を使えるようにしておきます。

切片があるパターン（切片なしを強制的に設定しない）で

lm.R2.RR(intercept = "existing") %>% plot()

で計算すると

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

となり、切片がありのパターンだと、R27がR21より概ね大きく、結構ばらついていることがわかる。

実際には、Rのlmではこのパターンだと、R21が計算出力さるわけで、R21で当てはまりが悪いとされるものがR27で高く評価されることがあるが、これは実際には出力されることはない。

切片がないパターン（強制的に原点を通す設定）をすると

lm.R2.RR(intercept = "zero") %>% plot()

で計算すると

f:id:indenkun:20200311171844p:plain — fig.2

で切片なしだと、R27とR21の差はfig.1と比較するとそれほど大きくないように見える。

こちらは実際にはR27が計算されるので、個人的な印象としてはそれほどR21と大きな差が出てこないようになっているのかも、と思った。

Kvalsethの式で計算したR²の差

計算用に使用した自作関数が、

lm.R2.kk <- function(n = 10, 
                     times = 100, 
                     intercept = c("existing", "zero")){
  kavl_R21 <- NULL
  kavl_R27 <- NULL
  lmintercept <- intercept
  for(i in 1:times){
    x <- rnorm(n)
    y <- rnorm(n)
    kavl_R21[i] <- lm.kavl_R2(x,y,R2 = "R21", intercept = lmintercept)
    kavl_R27[i] <- lm.kavl_R2(x,y,R2 = "R27", intercept = lmintercept)
  }
  data.frame(kavl_R21, kavl_R27)
}

lm用のものと設定は同じ。

切片ありのパターンで計算すると、

lm.R2.kk(intercept = "existing") %>% plot()

で

f:id:indenkun:20200311172805p:plain — fig.3

となり、fig.1と似たような結果になる。

R27のほうが適合度がいいような感じに出力されるが、一般的な決定係数R²の定義式がR21なので、当てはまりがよく評価されすぎることは少ない（？）のかもしれない。

切片なしのパターンで計算すると、

lm.R2.kk(intercept = "zero") %>% plot()

で

f:id:indenkun:20200311173140p:plain — fig.4

で、R21が負の値をとり得るので、概ねR27とR21は似たような値を取りつつ、おそらく切片を0にするのが非常に妥当でない場合にR21が負の値をとってひどく両者の値がずれていると思う。

R21が負の値を取らなければ一見fig.2と似ているので、これをもって切片がないときはR21だと負の値を取るからR27で計算しているよ（疑似的にだけど）、ということをlmでやっているのかもしれないと思った。

Kvalsethとlmの式でのRの差

Kvalsethとlmで同じR21、またはR27がどうなっているのかは次の自作関数で計算した。

lm.R2.RK <- function(n = 10,
                     times = 100,
                     RR2 = c("R21","R27"),
                     kavlR2 = c("R21","R27"),
                     intercept = c("existing","zero"))
{
  R_R2 <- NULL
  kavl_R2 <- NULL
  if(intercept == "existing"){
    lmform <- y~x
    lmintercept <- "existing"
  }
  else{
    lmform <- y~x-1
    lmintercept <- "zero"
  }
  for(i in 1:times){
    x <- rnorm(n)
    y <- rnorm(n)
    R_R2[i] <- lm.R_R2(form = lmform,R2 = RR2)
    kavl_R2[i] <- lm.kavl_R2(x,y, R2 = kavlR2, intercept = lmintercept)
  }
  if(RR2 == "R21"){
    R_R21 <- R_R2
    if(kavlR2 == "R21"){
      kavl_R21 <- kavl_R2
      dat <- data.frame(R_R21,kavl_R21)
    }
    else if(kavlR2 == "R27"){
      kavl_R27 <- kavl_R2
      dat <- data.frame(R_R21,kavl_R27)
    }
  }
  if(RR2 =="R27"){
    R_R27 <- R_R2
    if(kavlR2 == "R21"){
      kavl_R21 <- kavl_R2
      dat <- data.frame(R_R27,kavl_R21)
    }
    else if(kavlR2 == "R27"){
      kavl_R27 <- kavl_R2
      dat <-data.frame(R_R27,kavl_R27)
    }
  }
  dat
}

いままでのものに、RR2=でRの出力方式でRのR21とR27のいずれを計算するか、kavlR2=でR21とR27のいずれかを計算するかを指定できるようにした。

R21同士を切片ありで計算すると

lm.R2.RK(RR2 = "R21", kavlR2 = "R21", intercept = "existing") %>% plot()

で、

f:id:indenkun:20200311174131p:plain — fig.5

となり、RのlmはKvalsethのR21で求めてるよといっているので（実際の挙動はちょっと違う）、同じものを求めているのでR_R21も、kavl_R21も同じものが出る（当たり前の安心感）。

ちょっとずれているところがあるのは、浮動小数点の処理や丸めの処理のズレ程度？

R21同士を切片なしで計算すると

lm.R2.RK(RR2 = "R21", kavlR2 = "R21", intercept = "zero") %>% plot()

f:id:indenkun:20200311174713p:plain — fig.6

となり、試しに10000回施行とすると

lm.R2.RK(times = 10000, RR2 = "R21", kavlR2 = "R21", intercept = "zero") %>% plot()

f:id:indenkun:20200311174902p:plain — fig.7

となり、同じものを求めているとは思えない結果になる。

kavl_R21が負の値を取るのと、R_R21が正の値しか取らない関係での差が出ている。

概ね、R_R21で当てはまりが悪いとされる0に近いほど、kavl_R21でも大きく負の値をとっている傾向にあるが、そうじゃないものもるので一概には言い切れない。

R27同士を切片なしで計算すると

lm.R2.RK(RR2 = "R27", kavlR2 = "R27", intercept = "existing") %>% plot()

f:id:indenkun:20200311175536p:plain — fig.8

でkavl_R27もR_R27も同じ値をとっている。

R27同士で切片ありで計算すると

lm.R2.RK(RR2 = "R27", kavlR2 = "R27", intercept = "zero") %>% plot()

f:id:indenkun:20200311175942p:plain — fig.9

で、lmのD documentで同じ計算式だよといっているだけあって（実際には挙動は異なる）、結果として同じ物が出る（よかった）。

結局何を言いたいのか

原点を無理やり通したもののR²を計算するときは要注意（当たり前）。

そんな当たり前のことはおいておていて、実際のKvalsethの式とRのlmで計算されるR² は切片の有無で挙動が変わることでその出力される結果が一致することとになっていることがわかった。

実務的には問題はそう多くないと思う。

一方で、RがKvalsethの式と計算の挙動が異なることで、同じ決定係数R²を求めているのにこんなにも大きな差が出る可能性があることを考慮しておくのも一考だと思った。

備忘ログ

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

Rのlmの決定係数R^2について定義式の違いでの差がどの程度なのかプロットしてみた

これまでの簡単なまとめ

Rのlmで求める式でのR21とR27の差について

Kvalsethの式で計算したR²の差

Kvalsethとlmの式でのRの差

結局何を言いたいのか

これまでの簡単なまとめ

Rのlmで求める式でのR21とR27の差について

Kvalsethの式で計算したR2の差

Kvalsethとlmの式でのRの差

結局何を言いたいのか

Kvalsethの式で計算したR²の差