この記事はバージニア大学のhttps://data.library.virginia.edu/is-r-squared-useless/を抄訳したものである。
決定係数は役立たずなのか？
2015年10月、「講義で教授が『決定係数は役立たず』と言い放っていたがマジか」という投稿がRedditで注目を集めた。
この教授はカーネギーメロン大学のCosma Shalizi先生だとわかり、問題の講義の内容も、http://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/10/lecture-10.pdfの3.2節から書かれていると特定できた。
決定係数は、回帰モデルが目的変数の(平均からの)変動をどれだけ説明できているかを表す統計量である。決定係数が0.65なら、説明変数によって目的変数の変動の65%を説明できているということになる。
この理屈でいけば、決定係数の高いモデルが好ましいことになるが、Shalizi先生はこれに異議を唱えた。もちろん確たる根拠があってのことだ。
それでは、Shalizi先生が決定係数について主張したことをいくつか確認し、それぞれをRで書いたシミュレーションで検証してみよう。
【主張1】
決定係数は、モデルの当てはまりの良さを測るものではない。
モデルが完全に正しいときでも、その決定係数は任意に低くなりうるからだ。
例えば、単回帰モデルのすべての前提があらゆる点で正しく成り立っていても、その残差の従う正規分布の分散σ^2が大きくなるほど、決定係数は0に近づく1。
この主張を実証するのは簡単で、下記のRのコードを実行して確認できる。
ここでは、(1)線形回帰が成り立つ仮定を満たすデータを生成し、(2)データに単回帰を当てはめ、(3)決定係数を報告する「r2.0」関数を作成している。ただし簡便のため、関数引数は誤差の標準偏差σのみとしている。また(1)で作ったデータとは、説明変数xの個々の値に同じ1次関数をあてはめ作った目的変数yの値に、分散一定の正規分布からランダムに生成した誤差を足しあわせたものである。
コードでは、σを小さい順に並べた配列の個々の要素に「r2.0」関数を適用(apply)し、その結果をプロットしている。
r2.0 <- function(sig){
x <- seq(1,10,length.out = 100) # our predictor
y <- 2 + 1.2*x + rnorm(100,0,sd = sig) # our response; a function of x plus some random noise
summary(lm(y ~ x))$r.squared # print the R-squared value
}
sigmas <- seq(0.5,20,length.out = 20)
rout <- sapply(sigmas, r2.0) # apply our function to a series of sigma values
plot(rout ~ sigmas, type="b")
プロットしたグラフを見ると、線形モデルの仮定があらゆる面で完全に成り立つ状況にも関わらず、σが大きくなると決定係数が大きく減少していることがわかる。
【主張2】
モデルの前提が全くもって成り立っていなくても、決定係数は1に近い値を取ることがある。
これも、決定係数ではモデルの当てはまりの良さを測れない理由を示したものである。
Shalizi先生の講義ノートにある下記のコードを使って、説明変数xと目的変数yの関係が非線形になっているデータを作り、それを線形回帰してみよう。
set.seed(1)
x <- rexp(50,rate=0.005) # our predictor is data from an exponential distribution
y <- (x-1)^2 * runif(50, min=0.8, max=1.2) # non-linear data generation
plot(x,y) # clearly non-linear
xとyをプロットしたグラフは以下のようになる。
ここで決定係数を求めてみる。
summary(lm(y ~ x))$r.squared
# [1] 0.8485146
すると「xとyが線形関係にあること」というモデルの前提が全く成り立っていないのに、決定係数は約0.85という非常に高い値を取ってしまった。
この場合、決定係数の高さをもって良いモデルだと言い張ることはできない。
望ましくは先にデータをプロットし、このケースに単回帰を当てはめるのが不適切であると気づくべきであった。
【主張3】
決定係数は予測誤差について何も説明できない。
求まる係数と残差の分散がいずれも等しくなるデータを回帰分析したとき、その平均二乗誤差は常に等しくなる。しかし決定係数は、説明変数の値のとる範囲が変わると、0から1までの任意の値を取りえてしまう。
予測誤差を測定したいなら平均二乗誤差を使ったほうが良い。
下記のRのコードで説明変数xは、1から10まで等間隔に取られた100個のデータになっている。これを動かしてできた回帰式による決定係数は約0.94、平均二乗誤差は約0.65となる。
x <- seq(1,10,length.out = 100)
set.seed(1)
y <- 2 + 1.2*x + rnorm(100,0,sd = 0.9)
mod1 <- lm(y ~ x)
summary(mod1)$r.squared
# [1] 0.9383379
sum((fitted(mod1) - y)^2)/100 # Mean squared error
# [1] 0.6468052
他方、説明変数xの範囲を1から2までに狭めた以外は上とまったく同じコードを実行してみる。決定係数は約0.15、平均二乗誤差は約0.65となる。
x <- seq(1,2,length.out = 100)
set.seed(1)
y <- 2 + 1.2*x + rnorm(100,0,sd = 0.9)
mod1 <- lm(y ~ x)
summary(mod1)$r.squared
# [1] 0.1502448
sum((fitted(mod1) - y)^2)/100 # Mean squared error
# [1] 0.6468052
求まった平均二乗誤差は前と同じままだが、決定係数は約0.94から0.15まで下がってしまった。
求めた2つの回帰式の平均二乗誤差が変わらないということは、両者の予測能力は同じということになる。
しかし決定係数を見ると、値の高いモデルの方がより高い予測能力を持っているかのようにミスリードされてしまう。
【主張4】
目的変数yをそのまま学習に使ったモデルと、対数変換などを施したyを学習に使ったモデルを比較するのに決定係数は使えない。変換方法の相異なるyを使ったモデル同士についても同様だ。
また、モデルの前提がより良く満たされるようにyを変換した場合でさえ、あっさり決定係数が下がってしまうことがある。
線形回帰の前提を満たすために対数変換が必要となるデータを使って、このことを検証してみよう。
以下のRコードは先のものと違って、yの値をさらに指数関数化している。
これを実行したとき、回帰式の決定係数は約0.003と非常に低くなる。また残差プロットでは外れ値がいくつか確認され、線形回帰の前提のひとつである等分散性2が満たせていないことがわかる。
x <- seq(1,2,length.out = 100)
set.seed(1)
y <- exp(-2 - 0.09*x + rnorm(100,0,sd = 2.5))
summary(lm(y ~ x))$r.squared
# [1] 0.003281718
plot(lm(y ~ x), which=3)
この問題の解決にはデータを対数変換するのが一般的なため、それを下記のコードで試す。すると、今度は等分散性がより満たされた残差プロットが得られる。
しかし決定係数は約0.0007となり、以前より更に1桁下がってしまう！
plot(lm(log(y)~x),which = 3)
summary(lm(log(y)~x))$r.squared
# [1] 0.0006921086
これは極端な例で、いつもこのようになるとは限らない。実際、対数変換を行うと、通常は決定係数が増加する。
しかし先ほど示したように、線形回帰の前提をよりよく満たすようにしても、必ずしも高い決定係数が得られるとは限らない。
よって、決定係数はモデル間で比較できないといえる。
【主張5】
決定係数は、「回帰モデルが目的変数の分散をどれだけ説明できているかの割合」と言われることが非常に一般的である。
しかし単回帰の場合は、逆に目的変数で説明変数を予測する回帰式を作っても、同じ決定係数が得られる。
つまり高い決定係数が得られたとしても、ある変数を別の変数でどれだけ説明できるかについては何ともいえないのだ。
このことは下記のコードで簡単に確認できる。どちらも同じ約0.71という決定係数を取る。
x <- seq(1,10,length.out = 100)
y <- 2 + 1.2*x + rnorm(100,0,sd = 2)
summary(lm(y ~ x))$r.squared
# [1] 0.7065779
summary(lm(x ~ y))$r.squared
# [1] 0.7065779
この結果を見ると、xがyを説明しているのか、それとも逆にyがxを説明しているのかわからなくなってくる。我々は「原因」という言葉の周りでぐるぐる踊るために、「説明」という言葉を発しているのだろうか？とさえ思えてしまう。
しかし、このような単回帰分析において、決定係数は単なるxとyの相関の2乗にすぎない。そのことは、下記のコードで確認できる。
all.equal(cor(x,y)^2, summary(lm(x ~ y))$r.squared, summary(lm(y ~ x))$r.squared)
# [1] TRUE
ならば、決定係数の代わりに相関係数を使ってみてはどうか？しかし、相関係数も線形関係を要約するものに過ぎず、あらゆるデータの説明に適しているとは限らない。
この場合もやはり、データをプロットして説明変数と目的変数がどのような関係になっていそうか確認することが強く推奨される。
まとめ
これまでの話を要約しよう。
- 決定係数ではモデルの当てはまりの良さを測れない。
- 決定係数では予測誤差も測れない。(注: 測りたいなら平均二乗誤差などを使おう)
- 目的変数を変換してから学習したモデルと、変換せず学習したモデルとを、決定係数では比較できない。
- 決定係数は、ある変数が他の変数をどのように説明するかを測るものではない。
これまで取り上げた以外にも、Shalizi先生の講義ノートでは決定係数が無意味である理由をさらに多く述べている。
また調整済み決定係数でも、上記の問題のどれにも対処できないことも付記しておく必要があろう。
ならば、決定係数を使う理由はあるのだろうか？
Shalizi先生は「ない。決定係数が間違いなく役立つといえる状況を見たことがない」と断言している。
もちろん、統計学者やReddit読者の中にはそう思わない人もいるだろう。しかし、あなたの見解がどうであれ、データ分析の結果を説明するのに決定係数を使おうとするなら、あなたが伝えたいと思っていることがそれで実際に伝わるか、改めて検証するのが賢明だろう。