回帰分析の結果はどうなってるか

Google Earth Engineで使える線形回帰は４種類

GEEのガイドでは次の4種類のメソッドを使った線形回帰が紹介されている。これまでも使ったけど、係数以外の返り値がどうなってるのか詳しく見たい、似たように見えるけど、どこが違うのか見たいというのがやりたいこと。

ee.Reducer.linearFit()
ee.Reducer.linearRegression()
ee.Reducer.robustLinearRegression()
ee.Reducer.ridgeRegression()

簡単な違いの説明

単回帰分析 ee.Reducer.linearFit()

単回帰分析専門。引数はX（説明変数のバンド）とY（目的変数バンド）を指定すればOK。返り値は'scale'（傾き）と'offset'（切片）の２つのバンドを持つImage。

重回帰分析もできる ee.Reducer.linearRegression()

一般線形モデル。引数はX（説明変数のバンド数）とY（目的変数のバンド数）の二つ。デフォルトでXに定数項が入っていないので含めるのを忘れないこと。返り値は2バンドのArray Imageで、1バンド目の ['coefficients'] に推定された切片と係数が、2バンド目の ['residuals'] に目的変数のRMSEが入っている。

決定係数や検定の結果は入っていない。決定係数はRMSEと分散から計算できる↓

// R squared from RMSE
var rSquared = function(residualImage, coefficientsImage, imageCollection){
  var filteredIC = imageCollection.select('celsius');

  var n = filteredIC.count();
  var p = ee.Image.constant(coefficientsImage.bandNames().length().subtract(1));
  var one = ee.Image.constant(1);

  var variance = filteredIC.reduce(ee.Reducer.variance()); // .sampleVariance()だと不偏標本分散のようで値が異なる
  var rmseSq = residualImage.multiply(residualImage);

  var rSquared = one.subtract(rmseSq.divide(variance)); // 1 - rmseの二乗/分散
  var adjRSquared = ((n.subtract(one)).multiply(rSquared).subtract(p)).divide (n.subtract(p).subtract(one)); // ややこしい回帰のときはこれじゃあ力不足かも。

  return residualImage.addBands(rSquared)
                        .addBands(adjRSquared)
                        .rename(['RMSE','R-sqared','adjusted R-squared']);
};

これでRのlmでのsummaryの出力と一致。決定係数と自由度調整済み決定係数の計算式はwikipediaと他のブログに教えてもらいました。

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

Rのsummary関数を使わずに、決定係数を求めてみる | while(isプログラマ)

Xに説明変数をn乗した項を含めて多項式回帰にしたり、その他三角関数や対数関数の回帰も可能。気温データを使ったここまでの練習コードはこちら https://code.earthengine.google.com/1c414c3696fa60bf1e620715f2ea8faf

ロバスト回帰分析 ee.Reducer.robustLinearRegression()

外れ値の影響を受けにくいロバスト回帰分析。ほとんどee.Reducer.linearRegression()と同じだけど、引数にベータが加わるらしい。詳しくは論文を読まないとわからなそう（読んでもわからないかも）だったので、調べていません。

リッジ回帰分析 ee.Reducer.ridgeRegression()

リッジ回帰の引数はX（説明変数のバンド数）、Y（目的変数のバンド数）、ラムダ（正則化項のパラメータ、デフォルト値0.1）の３つ。.linearRegression()と違い、Xに定数項を含める必要はない。返り値は3バンドのArray Imageで、1バンド目の ['coefficients'] に推定された切片と係数が、2バンド目の ['residuals'] にはRMSEではなく残差平方和の平方根¹が入っている。3バンド目の['pValue'] は偏回帰係数をt検定したp値²が入っている。

ラムダを0にすれば、普通の線形回帰と同じ。なので、t検定の結果が欲しい場合は.linearRegression()ではなく.ridgeRegression()を使えば出てくる。residulasの出力がよくわかんないけど、残差平方和の平方根からも分散を使えば決定係数は計算できる。だいたいは.linearRegression()と同じなので省略、こちらをどうぞ↓

ridgeRegression()でλを0にした場合：https://code.earthengine.google.com/6f256e27b04e5dcc252f8ab9363ac90d

その他の回帰的なもの

チュートリアルにノンパラメトリックな方法が紹介されている。ランダムフォレスト回帰とかできるのかもしれないけど、探せなかった。調べたらできるのかも。 developers.google.com

回帰の分散分析（F検定）はやりかたがわからない

調べてもわからなかったのが、回帰式自体が有意かどうかのF検定の方法。F値と自由度は求められるけど、その値がF分布でどうなのよという部分がわからない。自作するとしたらベータ分布やらガンマ分布やら出てきてちょっと…という感じ。先人の知恵があるのか。そのうち調べたい³。

リファレンスを読むと"Additional outputs are a vector of the root mean square of the residuals of each dependent variable …“と書いてあるけど、Rで確かめたところmeanにはなっていないみたいで謎。↩
数学的な言い回しがわかんないんだけど、係数のt検定ということを言いたい。↩
さっき見たら ee.Image.erf() というのがあったんだけど、これでいけるのかな？わかんない。↩