Add updated `xr_regression` function for multi-dimensional linear regression

[robbibt]

Proposed changes

This PR updates the older lag_linregress_3d function into a new and improved xr_regression function for calculating useful regression statistics between two multi-dimensional xarray datasets, including:

• Regression slope, standard error, intercept
• P-value
• Correlation and covariance
• R-squared

For example, the function can be used to calculate regressions between two 3D datasets (e.g. time, x, y), or between a 3D (time, x, y) dataset and a 1D dataset (time):

This PR includes tests that verify that the results produced by this function are identical to those produced by the scipy.stats.linregress function (including for "two-sided", "less" and "greater" alternative hypotheses).

Checklist

(Replace [ ] with [x] to check off)

• [ ] Notebook created using the DEA-notebooks template
• [x] Remove any unused Python packages from Load packages
• [x] Remove any unused/empty code cells
• [ ] Remove any guidance cells (e.g. General advice)
• [x] Ensure that all code cells follow the PEP8 standard for code. The jupyterlab_code_formatter tool can be used to format code cells to a consistent style: select each code cell, then click Edit and then one of the Apply X Formatter options (YAPF or Black are recommended).
• [x] Include relevant tags in the final notebook cell (refer to the DEA Tags Index, and re-use tags if possible)
• [x] Clear all outputs, run notebook from start to finish, and save the notebook in the state where all cells have been sequentially evaluated
• [ ] Test notebook on both the NCI and DEA Sandbox (flag if not working as part of PR and ask for help to solve if needed)
• [ ] If applicable, update the Notebook currently compatible with the NCI|DEA Sandbox environment only line below the notebook title to reflect the environments the notebook is compatible with
• [x] Check for any spelling mistakes using the DEA Sandbox's built-in spellchecker (double click on markdown cells then right-click on pink highlighted words). For example:

[cbur24]

@robbibt Awesome! I was just looking for something like this function. A possible enhancement: any interest in including options for robust regression instead of OLS? This can be especially important for satellite time-series regression where slopes can be influenced by outlier values. Theil-sen slopes with a Mann Kendall test is a good example of robust regression.

scipy theil-slopes: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mstats.theilslopes.html

An example wrapper implemented for xarray here: https://github.com/josuemtzmo/xarrayMannKendall/blob/master/xarrayMannKendall/xarrayMannKendall.py

An example from my work implementing the above wrapper, go down to In [11] https://github.com/cbur24/AusENDVI/blob/main/notebooks/analysis/Trends_in_Seasonality.ipynb

[robbibt]

@robbibt Awesome! I was just looking for something like this function. A possible enhancement: any interest in including options for robust regression instead of OLS? This can be especially important for satellite time-series regression where slopes can be influenced by outlier values. Theil-sen slopes with a Mann Kendall test is a good example of robust regression.

scipy theil-slopes: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mstats.theilslopes.html

An example wrapper implemented for xarray here: https://github.com/josuemtzmo/xarrayMannKendall/blob/master/xarrayMannKendall/xarrayMannKendall.py

An example from my work implementing the above wrapper, go down to In [11] https://github.com/cbur24/AusENDVI/blob/main/notebooks/analysis/Trends_in_Seasonality.ipynb

Thanks @cbur24! Expanding this to include multiple regression methods would be super neat - we could have a simple param like method= to change the internal func and still return a consistent set of outputs.

The current implementation is designed to be almost completely vectorised array maths, which makes it really quick and non-memory hungry. It would be neat to be able to also do that across different regression methods - at first glance the xarrayMannKendall implementation above seems pretty complex and maybe not something that could be easily switched to array maths. I wonder if a simple addition at this stage could be adding something like robust outlier detection, e.g. with MAD or RANSAC or similar? Not quite as good as a dedicated robust regression, but could make it a bit more useful for really noisy data...

[cbur24]

Totally agree that it would be complex addition, and from reading how seaborn/pandas does robust regression it seems as though most robust regression techniques are slow and/or memory intensive. In light of that, having the option for outlier detection and removal sounds like a great (lightweight) addition that would do 95 % of what robust regression does but without the overhead.

[robbibt]

@cbur24 Have added a simple MAD outlier detection function and params in xr_regression to apply it (https://github.com/GeoscienceAustralia/dea-notebooks/pull/1226/commits/753b3686e4d11909bbd0a52244409415171844d8) - I'm not sure it's the best approach as the outlier detection doesn't take account of the relationship between the two variables (it's univariate only). But it's there anyway as an experimental feature, so we can see if it's useful! 🙂 Definitely keen to expand this func to robust regression in the future though, so if you ever want to raise a PR, feel free!

[cbur24]

@robbibt I did some basic testing of this function today using a couple of janky non-datacube netcdfs. The function works well (and fast!), with the exception that the lag parameter was failing. The following code:

reg = xr_regression(ndvi, vpd, dim='time', lag_y=1)

results in this error:

ValueError: conflicting sizes for dimension 'latitude': length 1 on <this-array> and length 680 on {'longitude': 'longitude', 'latitude': 'latitude', 'time': 'time'} at line 946: cov = ((x - xmean) * (y - ymean)).sum(dim=dim) / (n)

Weirdly, when I run the same lagged function call but with dask, I don't receive an error but the result is an all-NaN array. I'm guessing it maybe has something to do with dimension alignment.

The input xarray datasets ndvi and vpd look like this:

[robbibt]

@robbibt I did some basic testing of this function today using a couple of janky non-datacube netcdfs. The function works well (and fast!), with the exception that the lag parameter was failing. The following code:

reg = xr_regression(ndvi, vpd, dim='time', lag_y=1)

results in this error:

ValueError: conflicting sizes for dimension 'latitude': length 1 on <this-array> and length 680 on {'longitude': 'longitude', 'latitude': 'latitude', 'time': 'time'} at line 946: cov = ((x - xmean) * (y - ymean)).sum(dim=dim) / (n)

Weirdly, when I run the same lagged function call but with dask, I don't receive an error but the result is an all-NaN array. I'm guessing it maybe has something to do with dimension alignment.

The input xarray datasets ndvi and vpd look like this:

Thanks @cbur24! I'll admit that the lag functionality was the one bit I didn't test in this re-write. 🙃 Good catch! I'll see if I can reproduce this on my end, otherwise I might grab the files you're using if they're sharable (will let you know). If the lag stuff proves too complicated, I'm tempted to just leave it out of the function and let users handle that themselves outside of the func.

[review-notebook-app[bot]]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[robbibt]

Thanks for the great feedback @cbur24! There's a few issues with the current code (most importantly, dropna in the lag code would drop every timestep with any nodata... even a single pixel!), but the more I think about this, the more I think this functionality doesn't belong in xr_regression... it's not super clear exactly how the dimensions are aligned after lagging, and I think there's a very high risk of accidently producing results that are invalid and yet really opaque and difficult to troubleshoot.

I think it's probably best that users apply lags externally to the function, and make sure their datasets are perfectly compatible and ready for comparison before passing them in for comparison. 🙂