線形回帰

Earth Engine には、リダクタを使用して線形回帰を実行する方法がいくつかあります。

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

最も単純な線形回帰レジューサーは linearFit() で、定数項を持つ 1 変数の線形関数の最小二乗推定値を計算します。線形モデリングに柔軟なアプローチを適用するには、独立変数と従属変数の数を変更できる線形回帰レジューサーのいずれかを使用します。linearRegression() は通常の最小二乗回帰（OLS）を実装します。robustLinearRegression() は、回帰残差に基づくコスト関数を使用して、データ内の外れ値の重みを反復的に減らします（O’Leary、1990）。ridgeRegression() は、L2 正則化を使用した線形回帰を行います。

これらのメソッドを使用した回帰分析は、ee.ImageCollection、ee.Image、ee.FeatureCollection、ee.List オブジェクトを削減する場合に適しています。次の例は、それぞれのアプリケーションを示しています。linearRegression()、robustLinearRegression()、ridgeRegression() はすべて同じ入力と出力の構造を持ちますが、linearFit() は 2 バンドの入力（X の後に Y）を想定しています。ridgeRegression() には追加のパラメータ（lambda、省略可）と出力（pValue）があります。

ee.ImageCollection

linearFit()

データは 2 バンドの入力画像として設定する必要があります。最初のバンドは独立変数で、2 番目のバンドは従属変数です。次の例は、気候モデルによって予測される将来の降水量の線形傾向（NEX-DCP30 データの 2006 年以降）の推定を示しています。従属変数は予測される降水量で、独立変数は時間です。linearFit() を呼び出す前に追加します。

コードエディタ（JavaScript）

// This function adds a time band to the image.
var createTimeBand = function(image) {
  // Scale milliseconds by a large constant to avoid very small slopes
  // in the linear regression output.
  return image.addBands(image.metadata('system:time_start').divide(1e18));
};

// Load the input image collection: projected climate data.
var collection = ee.ImageCollection('NASA/NEX-DCP30_ENSEMBLE_STATS')
  .filter(ee.Filter.eq('scenario', 'rcp85'))
  .filterDate(ee.Date('2006-01-01'), ee.Date('2050-01-01'))
  // Map the time band function over the collection.
  .map(createTimeBand);

// Reduce the collection with the linear fit reducer.
// Independent variable are followed by dependent variables.
var linearFit = collection.select(['system:time_start', 'pr_mean'])
  .reduce(ee.Reducer.linearFit());

// Display the results.
Map.setCenter(-100.11, 40.38, 5);
Map.addLayer(linearFit,
  {min: 0, max: [-0.9, 8e-5, 1], bands: ['scale', 'offset', 'scale']}, 'fit');

出力には、「オフセット」（切片）と「スケール」の 2 つの帯域が含まれています（このコンテキストでの「スケール」は線の傾斜を指し、多くのリデューサーに入力されるスケール パラメータ（空間スケール）とは異なります）。結果は、増加傾向の領域が青色、減少傾向の領域が赤色、傾向なしの領域が緑色で、図 1 のようになります。


図 1. 予測される降水量に適用される linearFit() の出力。降水量が増加すると予測される地域は青色で、降水量が減少すると予測される地域は赤色で示されます。

linearRegression()

たとえば、2 つの従属変数（降水量と最高気温）と 2 つの独立変数（定数と時間）があるとします。このコレクションは前の例と同じですが、減算の前に定数バンドを手動で追加する必要があります。入力の最初の 2 つのバンドは「X」（独立）変数で、次の 2 つのバンドは「Y」（従属）変数です。この例では、まず回帰係数を取得し、次に配列画像をフラット化して、対象のバンドを抽出します。

コードエディタ（JavaScript）

// This function adds a time band to the image.
var createTimeBand = function(image) {
  // Scale milliseconds by a large constant.
  return image.addBands(image.metadata('system:time_start').divide(1e18));
};

// This function adds a constant band to the image.
var createConstantBand = function(image) {
  return ee.Image(1).addBands(image);
};

// Load the input image collection: projected climate data.
var collection = ee.ImageCollection('NASA/NEX-DCP30_ENSEMBLE_STATS')
  .filterDate(ee.Date('2006-01-01'), ee.Date('2099-01-01'))
  .filter(ee.Filter.eq('scenario', 'rcp85'))
  // Map the functions over the collection, to get constant and time bands.
  .map(createTimeBand)
  .map(createConstantBand)
  // Select the predictors and the responses.
  .select(['constant', 'system:time_start', 'pr_mean', 'tasmax_mean']);

// Compute ordinary least squares regression coefficients.
var linearRegression = collection.reduce(
  ee.Reducer.linearRegression({
    numX: 2,
    numY: 2
}));

// Compute robust linear regression coefficients.
var robustLinearRegression = collection.reduce(
  ee.Reducer.robustLinearRegression({
    numX: 2,
    numY: 2
}));

// The results are array images that must be flattened for display.
// These lists label the information along each axis of the arrays.
var bandNames = [['constant', 'time'], // 0-axis variation.
                 ['precip', 'temp']]; // 1-axis variation.

// Flatten the array images to get multi-band images according to the labels.
var lrImage = linearRegression.select(['coefficients']).arrayFlatten(bandNames);
var rlrImage = robustLinearRegression.select(['coefficients']).arrayFlatten(bandNames);

// Display the OLS results.
Map.setCenter(-100.11, 40.38, 5);
Map.addLayer(lrImage,
  {min: 0, max: [-0.9, 8e-5, 1], bands: ['time_precip', 'constant_precip', 'time_precip']}, 'OLS');

// Compare the results at a specific point:
print('OLS estimates:', lrImage.reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: ee.Geometry.Point([-96.0, 41.0]),
  scale: 1000
}));

print('Robust estimates:', rlrImage.reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: ee.Geometry.Point([-96.0, 41.0]),
  scale: 1000
}));

結果を調べると、linearRegression() の出力が linearFit() レジューサーによって推定された係数と同じであることがわかります。ただし、linearRegression() の出力には、他の従属変数 tasmax_mean の係数も含まれています。ロバストな線形回帰係数は OLS 推定値とは異なります。この例では、特定の時点でのさまざまな回帰方法の係数を比較しています。

ee.Image

ee.Image オブジェクトのコンテキストでは、回帰リダクタを reduceRegion または reduceRegions とともに使用して、領域内のピクセルに線形回帰を実行できます。次の例は、任意のポリゴン内の Landsat バンド間の回帰係数を計算する方法を示しています。

linearFit()

配列データのグラフについて説明しているガイドのセクションには、Landsat 8 SWIR1 バンドと SWIR2 バンドの相関関係の散布図が表示されます。ここでは、この関係の線形回帰係数が計算されます。'offset'（y 切片）と 'scale'（傾斜）のプロパティを含む辞書が返されます。

コードエディタ（JavaScript）

// Define a rectangle geometry around San Francisco.
var sanFrancisco = ee.Geometry.Rectangle([-122.45, 37.74, -122.4, 37.8]);

// Import a Landsat 8 TOA image for this region.
var img = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318');

// Subset the SWIR1 and SWIR2 bands. In the regression reducer, independent
// variables come first followed by the dependent variables. In this case,
// B5 (SWIR1) is the independent variable and B6 (SWIR2) is the dependent
// variable.
var imgRegress = img.select(['B5', 'B6']);

// Calculate regression coefficients for the set of pixels intersecting the
// above defined region using reduceRegion with ee.Reducer.linearFit().
var linearFit = imgRegress.reduceRegion({
  reducer: ee.Reducer.linearFit(),
  geometry: sanFrancisco,
  scale: 30,
});

// Inspect the results.
print('OLS estimates:', linearFit);
print('y-intercept:', linearFit.get('offset'));
print('Slope:', linearFit.get('scale'));

linearRegression()

ここでは、前のセクションの linearFit と同じ分析が適用されますが、今回は ee.Reducer.linearRegression 関数を使用します。回帰画像は、定数画像と、同じ Landsat 8 画像の SWIR1 バンドと SWIR2 バンドを表す画像の 3 つの別々の画像から構成されます。任意のバンドセットを組み合わせて、ee.Reducer.linearRegression による領域縮小用の入力画像を作成できます。バンドは同じソース画像に属している必要はありません。

コードエディタ（JavaScript）

// Define a rectangle geometry around San Francisco.
var sanFrancisco = ee.Geometry.Rectangle([-122.45, 37.74, -122.4, 37.8]);

// Import a Landsat 8 TOA image for this region.
var img = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318');

// Create a new image that is the concatenation of three images: a constant,
// the SWIR1 band, and the SWIR2 band.
var constant = ee.Image(1);
var xVar = img.select('B5');
var yVar = img.select('B6');
var imgRegress = ee.Image.cat(constant, xVar, yVar);

// Calculate regression coefficients for the set of pixels intersecting the
// above defined region using reduceRegion. The numX parameter is set as 2
// because the constant and the SWIR1 bands are independent variables and they
// are the first two bands in the stack; numY is set as 1 because there is only
// one dependent variable (SWIR2) and it follows as band three in the stack.
var linearRegression = imgRegress.reduceRegion({
  reducer: ee.Reducer.linearRegression({
    numX: 2,
    numY: 1
  }),
  geometry: sanFrancisco,
  scale: 30,
});

// Convert the coefficients array to a list.
var coefList = ee.Array(linearRegression.get('coefficients')).toList();

// Extract the y-intercept and slope.
var b0 = ee.List(coefList.get(0)).get(0); // y-intercept
var b1 = ee.List(coefList.get(1)).get(0); // slope

// Extract the residuals.
var residuals = ee.Array(linearRegression.get('residuals')).toList().get(0);

// Inspect the results.
print('OLS estimates', linearRegression);
print('y-intercept:', b0);
print('Slope:', b1);
print('Residuals:', residuals);

'coefficients' プロパティと 'residuals' プロパティを含む辞書が返されます。'coefficients' プロパティは、ディメンションが（numX、numY）の配列です。各列には、対応する従属変数の係数が含まれます。この場合、配列は 2 行 1 列で、1 行目 1 列目は y 切片、2 行目 1 列目は傾斜です。'residuals' プロパティは、各従属変数の残差の二乗平均平方根のベクトルです。結果を配列としてキャストし、目的の要素をスライスするか、配列をリストに変換してインデックス位置で係数を選択して、係数を抽出します。

ee.FeatureCollection

Sentinel-2 と Landsat 8 SWIR1 の反射率の線形関係を知りたいとします。この例では、ポイントの特徴コレクションとしてフォーマットされたピクセルのランダム サンプルを使用して、関係を計算します。最小二乗法の最適化曲線とともに、ピクセルペアの散布図が生成されます（図 2）。

コードエディタ（JavaScript）

// Import a Sentinel-2 TOA image.
var s2ImgSwir1 = ee.Image('COPERNICUS/S2/20191022T185429_20191022T185427_T10SEH');

// Import a Landsat 8 TOA image from 12 days earlier than the S2 image.
var l8ImgSwir1 = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044033_20191010');

// Get the intersection between the two images - the area of interest (aoi).
var aoi = s2ImgSwir1.geometry().intersection(l8ImgSwir1.geometry());

// Get a set of 1000 random points from within the aoi. A feature collection
// is returned.
var sample = ee.FeatureCollection.randomPoints({
  region: aoi,
  points: 1000
});

// Combine the SWIR1 bands from each image into a single image.
var swir1Bands = s2ImgSwir1.select('B11')
  .addBands(l8ImgSwir1.select('B6'))
  .rename(['s2_swir1', 'l8_swir1']);

// Sample the SWIR1 bands using the sample point feature collection.
var imgSamp = swir1Bands.sampleRegions({
  collection: sample,
  scale: 30
})
// Add a constant property to each feature to be used as an independent variable.
.map(function(feature) {
  return feature.set('constant', 1);
});

// Compute linear regression coefficients. numX is 2 because
// there are two independent variables: 'constant' and 's2_swir1'. numY is 1
// because there is a single dependent variable: 'l8_swir1'. Cast the resulting
// object to an ee.Dictionary for easy access to the properties.
var linearRegression = ee.Dictionary(imgSamp.reduceColumns({
  reducer: ee.Reducer.linearRegression({
    numX: 2,
    numY: 1
  }),
  selectors: ['constant', 's2_swir1', 'l8_swir1']
}));

// Convert the coefficients array to a list.
var coefList = ee.Array(linearRegression.get('coefficients')).toList();

// Extract the y-intercept and slope.
var yInt = ee.List(coefList.get(0)).get(0); // y-intercept
var slope = ee.List(coefList.get(1)).get(0); // slope

// Gather the SWIR1 values from the point sample into a list of lists.
var props = ee.List(['s2_swir1', 'l8_swir1']);
var regressionVarsList = ee.List(imgSamp.reduceColumns({
  reducer: ee.Reducer.toList().repeat(props.size()),
  selectors: props
}).get('list'));

// Convert regression x and y variable lists to an array - used later as input
// to ui.Chart.array.values for generating a scatter plot.
var x = ee.Array(ee.List(regressionVarsList.get(0)));
var y1 = ee.Array(ee.List(regressionVarsList.get(1)));

// Apply the line function defined by the slope and y-intercept of the
// regression to the x variable list to create an array that will represent
// the regression line in the scatter plot.
var y2 = ee.Array(ee.List(regressionVarsList.get(0)).map(function(x) {
  var y = ee.Number(x).multiply(slope).add(yInt);
  return y;
}));

// Concatenate the y variables (Landsat 8 SWIR1 and predicted y) into an array
// for input to ui.Chart.array.values for plotting a scatter plot.
var yArr = ee.Array.cat([y1, y2], 1);

// Make a scatter plot of the two SWIR1 bands for the point sample and include
// the least squares line of best fit through the data.
print(ui.Chart.array.values({
  array: yArr,
  axis: 0,
  xLabels: x})
  .setChartType('ScatterChart')
  .setOptions({
    legend: {position: 'none'},
    hAxis: {'title': 'Sentinel-2 SWIR1'},
    vAxis: {'title': 'Landsat 8 SWIR1'},
    series: {
      0: {
        pointSize: 0.2,
        dataOpacity: 0.5,
      },
      1: {
        pointSize: 0,
        lineWidth: 2,
      }
    }
  })
);


図 2. Sentinel-2 と Landsat 8 の SWIR1 TOA 反射率を表すピクセルのサンプルの散布図と最小二乗線形回帰線。

ee.List

2 次元 ee.List オブジェクトの列は、回帰レジューサーへの入力として使用できます。次の例は簡単な証明を示しています。独立変数は、y 切片が 0 で傾斜が 1 となる従属変数のコピーです。

linearFit()

コードエディタ（JavaScript）

// Define a list of lists, where columns represent variables. The first column
// is the independent variable and the second is the dependent variable.
var listsVarColumns = ee.List([
  [1, 1],
  [2, 2],
  [3, 3],
  [4, 4],
  [5, 5]
]);

// Compute the least squares estimate of a linear function. Note that an
// object is returned; cast it as an ee.Dictionary to make accessing the
// coefficients easier.
var linearFit = ee.Dictionary(listsVarColumns.reduce(ee.Reducer.linearFit()));

// Inspect the result.
print(linearFit);
print('y-intercept:', linearFit.get('offset'));
print('Slope:', linearFit.get('scale'));

ee.Array に変換して転置し、ee.List に戻すことで、変数が行で表されている場合はリストを転置します。

コードエディタ（JavaScript）

// If variables in the list are arranged as rows, you'll need to transpose it.
// Define a list of lists where rows represent variables. The first row is the
// independent variable and the second is the dependent variable.
var listsVarRows = ee.List([
  [1, 2, 3, 4, 5],
  [1, 2, 3, 4, 5]
]);

// Cast the ee.List as an ee.Array, transpose it, and cast back to ee.List.
var listsVarColumns = ee.Array(listsVarRows).transpose().toList();

// Compute the least squares estimate of a linear function. Note that an
// object is returned; cast it as an ee.Dictionary to make accessing the
// coefficients easier.
var linearFit = ee.Dictionary(listsVarColumns.reduce(ee.Reducer.linearFit()));

// Inspect the result.
print(linearFit);
print('y-intercept:', linearFit.get('offset'));
print('Slope:', linearFit.get('scale'));

linearRegression()

ee.Reducer.linearRegression() の適用は、上記の linearFit() の例と似ていますが、定数独立変数が含まれています。

コードエディタ（JavaScript）

// Define a list of lists where columns represent variables. The first column
// represents a constant term, the second an independent variable, and the third
// a dependent variable.
var listsVarColumns = ee.List([
  [1, 1, 1],
  [1, 2, 2],
  [1, 3, 3],
  [1, 4, 4],
  [1, 5, 5]
]);

// Compute ordinary least squares regression coefficients. numX is 2 because
// there is one constant term and an additional independent variable. numY is 1
// because there is only a single dependent variable. Cast the resulting
// object to an ee.Dictionary for easy access to the properties.
var linearRegression = ee.Dictionary(
  listsVarColumns.reduce(ee.Reducer.linearRegression({
    numX: 2,
    numY: 1
})));

// Convert the coefficients array to a list.
var coefList = ee.Array(linearRegression.get('coefficients')).toList();

// Extract the y-intercept and slope.
var b0 = ee.List(coefList.get(0)).get(0); // y-intercept
var b1 = ee.List(coefList.get(1)).get(0); // slope

// Extract the residuals.
var residuals = ee.Array(linearRegression.get('residuals')).toList().get(0);

// Inspect the results.
print('OLS estimates', linearRegression);
print('y-intercept:', b0);
print('Slope:', b1);
print('Residuals:', residuals);