רגרסיה לינארית

ב-Earth Engine יש כמה שיטות לביצוע רגרסיה ליניארית באמצעות reducers:

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

הפונקציה הפשוטה ביותר להפחתת רגרסיה לינארית היא linearFit(), שמחשבת את האומדן של הריבועים הקטנים ביותר של פונקציה לינארית של משתנה אחד עם ביטוי קבוע. כדי לקבל גישה גמישה יותר לבניית מודלים לינאריים, אפשר להשתמש באחד מהמקטינים של רגרסיה ליניארית שמאפשרים מספר משתנים בלתי קבוע של משתנים עצמאיים ותלויים. linearRegression() מטמיע רגרסיה רגילה של הריבועים הקטנים ביותר(OLS). ב-robustLinearRegression() נעשה שימוש בפונקציית עלות שמבוססת על ערכים שיורית של רגרסיה כדי לבטל באופן איטרטיבי את המשקל של ערכים חריגים בנתונים (O'Leary, 1990). ridgeRegression() מבצעת רגרסיה לינארית עם רגולריזציה מסוג L2.

ניתוח רגרסיה בשיטות האלה מתאים לצמצום אובייקטים מסוג ee.ImageCollection,‏ ee.Image,‏ ee.FeatureCollection ו-ee.List. הדוגמאות הבאות מדגימות את השימוש בכל אחת מהן. שימו לב שלפונקציות linearRegression(), ‏ robustLinearRegression() ו-ridgeRegression() יש את אותן מבני קלט ופלט, אבל לפונקציה linearFit() נדרש קלט של שני פסים (X ואחריו Y), ולפונקציה ridgeRegression() יש פרמטר נוסף (lambda, אופציונלי) ופלט (pValue).

ee.ImageCollection

linearFit()

צריך להגדיר את הנתונים כתמונה קלט עם שני פסים, כאשר הפס הראשון הוא המשתנה הבלתי תלוי והפס השני הוא המשתנה התלוי. בדוגמה הבאה מוצגת הערכה של המגמה הליניארית של משקעים עתידיים (אחרי 2006 בנתוני NEX-DCP30) שחזתה על ידי מודלים אקלימיים. המשתנה התלוי הוא כמות המשקעים הצפויה, והמשתנה הבלתי תלוי הוא הזמן, שנוסף לפני הקריאה ל-linearFit():

Code Editor‏ (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');

שימו לב שהפלט מכיל שתי רצועות, 'offset' (נקודת החיתוך) ו-'scale' ('scale' בהקשר הזה מתייחס לשיפוע הקו, ואין להתבלבל בינו לבין פרמטר scale שמוזן למצמצמים רבים, שהוא ההיקף המרחבי). התוצאה, עם אזורים של מגמה עולה בכחול, מגמה יורדת באדום ואזורים ללא מגמה בירוק, אמורה להיראות בערך כמו באיור 1.


איור 1. הפלט של linearFit() שחלה על כמות המשקעים הצפויה. אזורים שבהם צפויה עלייה בכמות המשקעים מוצגים בכחול, ואזורים שבהם צפויה ירידה בכמות המשקעים מוצגים באדום.

linearRegression()

לדוגמה, נניח שיש שני משתנים תלויים: משקעים וטמפרטורה מקסימלית, ושני משתנים בלתי תלויים: קבוע וזמן. האוסף זהה לדוגמה הקודמת, אבל צריך להוסיף את הפס הקבוע באופן ידני לפני ההפחתה. שתי הפסקות הראשונות של הקלט הן המשתנים 'X' (בלתי תלויים), והשתיים הבאות הן המשתנים 'Y' (תלויים). בדוגמה הזו, קודם מקבלים את גורמי הרגרסיה, ואז משטחים את תמונת המערך כדי לחלץ את הפסקות העניין:

Code Editor‏ (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()

בקטע של המדריך שמתאר תרשימים של נתוני מערך מוצג תרשים פיזור של המתאם בין הפסגות SWIR1 ו-SWIR2 של Landsat 8. כאן מחושבים מקדמי הרגרסיה הליניארית של הקשר הזה. הפונקציה מחזירה מילון שמכיל את המאפיינים 'offset' (נקודת החיתוך עם ציר y) ו-'scale' (שיפוע).

Code Editor‏ (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. חשוב לזכור שתמונה של רגרסיה מורכבת משלושה תמונות נפרדות: תמונה קבועה ותמונות שמייצגות את הפס SWIR1 ואת הפס SWIR2 מאותה תמונה של Landsat 8. חשוב לזכור שאפשר לשלב כל קבוצה של פסים כדי ליצור תמונה של קלט לצורך צמצום אזור באמצעות ee.Reducer.linearRegression, והן לא חייבות להשתייך לאותה תמונה מקור.

Code Editor‏ (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). כל עמודה מכילה את המקדמים של המשתנה התלוי התואם. במקרה הזה, המערך מכיל שתי שורות ועמודה אחת: שורה אחת, עמודה אחת היא נקודת החיתוך עם ציר ה-y, ושורה שתיים, עמודה אחת היא השיפוע. המאפיין 'residuals' הוא הווקטור של השורש הריבועי הממוצע של השגיאות הנותרות של כל משתנה תלוי. כדי לחלץ את המקדמים, מעבירים את התוצאה למערך ואז מנתחים את הרכיבים הרצויים, או ממירים את המערך לרשימה ובוחרים את המקדמים לפי מיקום האינדקס.

ee.FeatureCollection

נניח שאתם רוצים לדעת מה הקשר הליניארי בין השתקפות SWIR1 של Sentinel-2 לבין השתקפות SWIR1 של Landsat 8. בדוגמה הזו, נעשה שימוש במדגם אקראי של פיקסלים בפורמט של אוסף תכונות של נקודות כדי לחשב את הקשר. נוצר גרף פיזור של זוגות הפיקסלים יחד עם קו המתאים ביותר של הכי פחות ריבועים (איור 2).

Code Editor‏ (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. תרשים פיזור וקו רגרסיה לינארי של הכי קטן בכיכר (least squares) לדגימת פיקסלים שמייצגים את ההחזרה של Sentinel-2 ו-Landsat 8 SWIR1 TOA.

ee.List

עמודות של אובייקטים ee.List דו-ממדיים יכולות לשמש כקלט למצמצמי רגרסיה. הדוגמאות הבאות מספקות הוכחות פשוטות. המשתנה הבלתי תלוי הוא עותק של המשתנה התלוי, שמניב ציר X חותך שווה ל-0 ושיפוע שווה ל-1.

linearFit()

Code Editor‏ (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.

Code Editor‏ (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()‎, מלבד העובדה שמשתנה בלתי תלוי קבוע נכלל.

Code Editor‏ (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);