Correlations and their significance with Pyleoclim

Correlations and their significance with `Pyleoclim`#

Authors#

Julien Emile-Geay, USC Earth Sciences

Feng Zhu, NUIST

Author1 = {“name”: “Julien Emile-Geay”, “affiliation”: “Department of Earth Sciences, University of Southern California”, “email”: “julieneg@usc.edu”, “orcid”: “0000-0001-5920-4751”}

Author2 = {“name”: “Feng Zhu”, “affiliation”: “Nanjing University of Information Science and Technology”, “email”: “fzhu@nuist.edu”, “orcid”: “0000-0002-9969-2953”}

Preamble#

Introduction#

Correlation analysis, despite its simplicity and many shortcomings, remains a centerpiece of empirical analysis in many fields, particularly the paleosciences. Computing correlations is trivial enough; the difficulty lies in properly assessing their significance. Of particular importance are three considerations:

Persistence, which violates the standard assumption that the data are independent (which underlies most classical test of significance implemented, e.g. in Excel, Matlab or basic scipy).
Time irregularities, for instance comparing two records with different time axes, possibly unevenly spaced (which standard software cannot deal with out of the box)
Test multiplicity aka the “Look Elsewhere effect”, which states that repeatedly performing the same test can result in unacceptably high type I error (accepting correlations as significant, when in fact they are not). This arises e.g. when correlating a paleoclimate record with an instrumental field, assessing significance at thousands of grid points at once, or assessing significance within an age ensemble.

Accordingly, Pyleoclim facilitates an assessment of correlations that deals with all these cases, makes the necessary data transformations transparent to the user, and allows for one-line plot commands to visualize the results.

Goals#

In this notebook, we illustrate some of these capabilities with several simple scenarios, including:

Correlating a single Series object with another, with different time axes
Correlating a single Series object with a gridded field (using MultipleSeries)
Correlating a single Series object with a EnsembleSeries
Correlating two EnsembleSeries objects together

For a primer on these data structures, see the documentation.

Reading Time: 20 min

Running Time: several hours due to large SST grid

Keywords#

Correlations, Surrogates, Significance, Visualization

Pre-requisites#

This tutorial assumes basic knowledge of Python and the Series and MultipleSeries objects in Pyleoclim. If you are not familiar with this coding language, check out this tutorial: http://linked.earth/ec_workshops_py/. some basics statistics would be helpful, but are not strictly necessary to proceed.

Relevant Packages#

pyleoclim, matplotlib, xarray, cartopy

Data Description and Processing#

Coral Sr/Ca series (a proxy for sea-surface temperature) by Nurhati et al (2011). The data are from Palmyra Island in the tropical Pacific, at the edge of the NINO3.4 box (a popular measure of the state of El Niño Southern Oscillation).
Sea-surface temperature analysis of ERSSTv5 (Huang et al, 2017)

import pyleoclim as pyleo
pyleo.set_style('web')
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np

First let’s load the coral data:

D = LiPD()

D.load('../data/Ocn-Palmyra.Nurhati.2011.lpd')

timeseries,df = D.get_timeseries(D.get_all_dataset_names(),to_dataframe=True)

df

Loading 1 LiPD files

100%|██████████| 1/1 [00:00<00:00, 63.50it/s]

Loaded..
Extracting timeseries from dataset: Ocn-Palmyra.Nurhati.2011 ...

	mode	time_id	geo_meanLon	geo_meanLat	geo_meanElev	geo_type	geo_siteName	geo_ocean	geo_pages2kRegion	hasUrl	...	paleoData_variableName	paleoData_hasMinValue	paleoData_variableType	paleoData_values	paleoData_dataType	paleoData_description	paleoData_inferredVariableType	paleoData_pages2kID	paleoData_ocean2kID	paleoData_inCompilation
0	paleoData	age	-162.13	5.87	-10.0	Feature	Palmyra	WP	Ocean	https://data.mint.isi.edu/files/lipd/Ocn-Palmy...	...	d18O	-5.75	measured	[-5.41, -5.47, -5.49, -5.43, -5.48, -5.53, -5....	NaN	NaN	NaN	NaN	NaN	NaN
1	paleoData	age	-162.13	5.87	-10.0	Feature	Palmyra	WP	Ocean	https://data.mint.isi.edu/files/lipd/Ocn-Palmy...	...	year	1886.13	inferred	[1998.29, 1998.21, 1998.13, 1998.04, 1997.96, ...	float	Year AD	Year	NaN	NaN	NaN
2	paleoData	age	-162.13	5.87	-10.0	Feature	Palmyra	WP	Ocean	https://data.mint.isi.edu/files/lipd/Ocn-Palmy...	...	Sr_Ca	8.78	measured	[8.96, 8.9, 8.91, 8.94, 8.92, 8.89, 8.87, 8.81...	NaN	NaN	NaN	Ocn_129	PacificNurhati2011	Ocean2k_v1.0.0
3	paleoData	age	-162.13	5.87	-10.0	Feature	Palmyra	WP	Ocean	https://data.mint.isi.edu/files/lipd/Ocn-Palmy...	...	d18O	0.07	measured	[0.39, 0.35, 0.35, 0.35, 0.36, 0.22, 0.33, 0.3...	NaN	NaN	NaN	NaN	NaN	NaN
4	paleoData	age	-162.13	5.87	-10.0	Feature	Palmyra	WP	Ocean	https://data.mint.isi.edu/files/lipd/Ocn-Palmy...	...	year	1886.13	inferred	[1998.21, 1998.13, 1998.04, 1997.96, 1997.88, ...	float	Year AD	Year	NaN	NaN	NaN

5 rows × 93 columns

row = df[df['paleoData_variableName']=='Sr_Ca']

srca = pyleo.Series(
    time=row['year'].values[0],
    value=row['paleoData_values'].values[0],
    time_name = 'Years',
    time_unit = row['yearUnits'],
    value_name = 'Sr_Ca',
    value_unit = row['paleoData_units'],
    label='Ocn-Palmyra.Nurhati.2011'
) # load only the Sr/Ca observations
srca.plot()

Time axis values sorted in ascending order

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Years [2    AD\nName: yearUnits, dtype: object]', ylabel='Sr_Ca [2    mmol/mol\nName: paleoData_units, dtype: object]'>)

../_images/e13845ad6a73a0bf55737728fcc5d52fb8a647333cc483ba8a545f81d21c2433.png

Next we load the ERSSTv5 data. Note that the actual dat are not being transferred over the network yet; all we are getting is a pointer towards a netCDF file on a NOAA THREDDS server. The real crunching will come later.

Let’s select a smaller portion of this dataset to correlate with the Sr/Ca record:

And let’s plot that domain:

fig = plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.Robinson(central_longitude=180))
ax.coastlines()
ax.gridlines()
ersstv5_tp.sst.isel(time=-1).plot(
    ax=ax, transform=ccrs.PlateCarree(), vmin=16, vmax=30, cmap='coolwarm'
)

<cartopy.mpl.geocollection.GeoQuadMesh at 0x16896f4c0>

../_images/65ea99ab1b1e73b64e2c434f1015f046e6e0f154ad0ae80b420fef3aa22e5bba.png

Now extract the NINO3.4 index:

nino34_box = ersstv5_tp.sel(lat=slice(5, -5), lon=slice(190, 240))
gb = nino34_box.sst.groupby('time.month')
nino34_anom = gb - gb.mean(dim='time')
weights = np.cos(np.deg2rad(nino34_box.lat))
nino34 = nino34_anom.weighted(weights).mean(dim=['lat', 'lon'])
nino34.plot()

[<matplotlib.lines.Line2D at 0x168ee9b80>]

../_images/968e9c43f6edb298cd4f437aa4526598529fd61ae02fa417413e734fb6dcae69.png

The date format is datetime64, which is not understood by Pyleoclim. We thus need to transform it to decimal years, thanks to this little trick.

dates = nino34['time'].values
years = dates.astype('datetime64[Y]').astype(int) + 1970
months = dates.astype('datetime64[M]').astype(int) % 12 + 1
decimal_year = years + months/12.0

nino34ts = pyleo.Series(time = decimal_year, value = nino34.values, value_name='NINO3.4',
                        value_unit= '$^\circ$C', label='Niño3.4 SST')
nino34ts.plot()

(<Figure size 720x288 with 1 Axes>,
 <AxesSubplot:xlabel='time', ylabel='NINO3.4 [$^\\circ$C]'>)

../_images/697e2a8e5efdbdfd9e5168c1faeae570c8f676405e5fbf6a46a4065886536bc3.png

Now we can compute correlations between the Nurhati et al (2011) Sr/Ca dataset and NINO3.4.

Correlations with `EnsembleSeries` objects#

Much like correlations between a Series and a MultipleSeries object, one can correlate Series and EnsembleSeries objects. To quote the documentation: The EnsembleSeries object is a child of the MultipleSeries object, that is, a special case of MultipleSeries, aiming for ensembles of related series. Ensembles usually arise from age modeling or Bayesian calibrations. All members of an EnsembleSeries object are assumed to share identical labels and units.

In our field correlation case, we could have treated the grid points as an EnsembleSeries object, since they share the same time axis and units. However, we prefer to illustrate EnsembleSeries on a case more likely to arise in paleoclimatology: an ensemble of realizations from a noise model.

Let’s use the gen_ts function to generate some fractal noise (this is the default setting) and add it to the Sr/Ca record:

n_ens = 100  # ensemble size
ts_list = []
for i in range(n_ens):
    _, noise = pyleo.utils.gen_ts(t=srca.time, model= 'colored_noise') 
    z, _, _ = pyleo.utils.standardize(noise)
    new_ts = pyleo.Series(time=srca.time, value=srca.value+z/10)
    ts_list.append(new_ts)
    
srcaEns = pyleo.EnsembleSeries(ts_list) # create ensemble series object

EnsembleSeries objects come with two methods unavailable to MultipleSeries objects: plot_traces() and plot_envelope(), the latter visualizing the quantiles of the ensemble through time.

fig, ax = srcaEns.plot_envelope()
srca.plot(color='k',ax=ax, zorder=100, lgd_kwargs={'bbox_to_anchor': (1, 1)}) # plot the original series on top

<AxesSubplot:xlabel='Year [AD]', ylabel='Sr_Ca [mmol/mol]'>

../_images/f23a270002e7fae92dc41575ff95d2eb138a15e90dc111736faf7804d115eac7.png

corrEns = srcaEns.correlation(srca)
corrEns.plot(title='Correlation Ensemble')

Looping over 100 Series in the ensemble

100%|█████████████████████████████████████████| 100/100 [00:43<00:00,  2.28it/s]

(<Figure size 288x288 with 1 Axes>,
 <AxesSubplot:title={'center':'Correlation Ensemble'}, xlabel='$r$', ylabel='Count'>)

../_images/e4e677f8b1d5ff5efe7aee3914b1100a8ff1c0b7f1cc620bb16fed37f1e8a6a2.png

Non-surprisingly, all correlations are positive and significant, by construction.

As with MultipleSeries, correlation() would also work on the EnsembleSeries object by itself (i.e. without specifying the target).

`EnsembleSeries` vs `EnsembleSeries`#

FInally, we tackle the case of correlating two ensemble series objects. As before, we turn NINO3.4 into an ensemble by adding some noise to it, and to spice it up we use a different ensemble size.

Suppose we name the two groups A and B, then the calculation of correlations will follow a pairwise procedure, i.e., compare the 1st (2nd, 3rd, …) member from A to the 1st (2nd, 3rd, ..) member from B, and stop when running out of members from either group.

n_ens = 50
ts_list = []
for i in range(n_ens):
    _, noise = pyleo.utils.gen_ts(t=nino34ts.time, model= 'colored_noise') 
    z, _, _ = pyleo.utils.standardize(noise)
    new_ts = pyleo.Series(time=nino34ts.time, value=nino34ts.value+z/10)
    ts_list.append(new_ts)
    
nino34Ens = pyleo.EnsembleSeries(ts_list) # create ensemble series object

corrEns2 = srcaEns.correlation(nino34Ens)
corrEns2.plot(title='Correlation: NINO3.4 vs Palmyra Sr/Ca Ensembles')

Looping over 100 Series in the ensemble

100%|█████████████████████████████████████████| 100/100 [00:44<00:00,  2.25it/s]

(<Figure size 288x288 with 1 Axes>,
 <AxesSubplot:title={'center':'Correlation: NINO3.4 vs Palmyra Sr/Ca Ensembles'}, xlabel='$r$', ylabel='Count'>)

../_images/7e51aea3f2fd2a6a2b3c89320b991c62b24ee4ce73b43931ce7d03ab77f02cd8.png

Takeaways#

pyleoclim allows to easily compute correlations between series with different time axes (including different resolutions).
pyleoclim assesses significance with 3 different methods that are meant to guard against the deleterious effects of persistence, thereby lowering the risk of spurious correlations and p-hacking.
pyleoclim allows to compute correlations between a series and a collection of series, e.g. coming from a gridded field, addressing the Multiple Hypothesis Problem with the False Discovery Rate.
pyleoclim allows to compute correlations between, and within, MultipleSeries and EnsembleSeries object, using those same tools.

Acknowledgments#

Project Pythia for this notebook on xarray manipulations

About `apply_fdr` in `EnsembleSeries.correlation()`#

One major difference between MultipleSeries.correlation() and EnsembleSeries.correlation() is that the latter has the argument apply_fdr=True by default, which means that in an ensemble context, the FDR approach is to applied by default to determine whether the calculated p-values are significant. You can opt out of it by setting the flag to False, but not before carefully considering whether this is wise.

Below we explain the basic idea of FDR and demonstrate the use of Pyleoclim’s underlying FDR function with two simple tests to help understand the idea.

Basic idea of FDR#

Quoting from Ventura et al. (2004): “A single test performed at signiﬁcance level \(\alpha\) has probability \(\alpha\) of rejecting the null hypothesis when it is in fact true. Hence if \(n\) such tests are performed when all \(n\) null hypotheses are true (the collective null hypothesis), then the average number of tests for which the null is falsely rejected is \(n\alpha\). For example, with \(\alpha = 5\%\), testing for a trend at 1000 locations at which no change really occurred would yield 50 signiﬁcant locations on average; this is unacceptably high.”

To overcome this disadvantage, there have been some alternatives, among which the enormously popular False Discovery Rate (FDR), introduced by Benjamini & Hochberg (1995) (BH95 hereafter). There are other variants such as Yekutieli & Benjamini (1999, 2001), and Ventura et al. (2004).

The idea is as follows: a multiple testing procedure, returns a list of \(p\)-values. With the traditional false positive test, the significance level \(\alpha=0.05\) is refered as the false positive rate (FPR), and we reject null hypothesis \(H_0\) at all locations \(p_i < \alpha\).

In the original BH95 method, the nominal FDR, denoted by \(q\), is defined as the rate we are willing to allow of false rejections out of all rejections, and is usually set to \(5\%\). The BH95 procedure rejects \(H_0\) in all cases \(i\) for which \(p_i \le p_k\), where \(k = \max_{i=0,...,n} \{i: p_{(i)}\le q i/n \}\). For an illustration, see Fig 2 of Hu et al (2017).

The underlying FDR function in `Pyleoclim`#

from pyleoclim.utils.correlation import fdr
?fdr

Signature: fdr(pvals, qlevel=0.05, method='original', adj_method=None, adj_args={})
Docstring:
Determine significance based on the false discovery rate

The false discovery rate is a method of conceptualizing the rate of type I errors in null hypothesis testing when conducting multiple comparisons. 
Translated from fdr.R by Dr. Chris Paciorek 

Parameters
----------

pvals : list or array
    A vector of p-values on which to conduct the multiple testing.

qlevel : float
    The proportion of false positives desired.

method : {'original', 'general'}
    Method for performing the testing.
        - 'original' follows Benjamini & Hochberg (1995);
        - 'general' is much more conservative, requiring no assumptions on the p-values (see Benjamini & Yekutieli (2001)).
        'original' is recommended, and if desired, using 'adj_method="mean"' to increase power.

adj_method: {'mean', 'storey', 'two-stage'}
    Method for increasing the power of the procedure by estimating the proportion of alternative p-values.
        - 'mean', the modified Storey estimator in Ventura et al. (2004)
        - 'storey', the method of Storey (2002)
        - 'two-stage', the iterative approach of Benjamini et al. (2001)

adj_args : dict
    Arguments for adj_method; see prop_alt() for description,
    but note that for "two-stage", qlevel and fdr_method are taken from the qlevel and method arguments for fdr()

Returns
-------

fdr_res : array or None
    A vector of the indices of the significant tests; None if no significant tests

See also
--------

pyleoclim.utils.correlation.corr_sig : Estimates the Pearson's correlation and associated significance between two non IID time series


References
----------

- fdr.R by Dr. Chris Paciorek: https://www.stat.berkeley.edu/~paciorek/research/code/code.html
File:      ~/opt/anaconda3/envs/paleodev/lib/python3.9/site-packages/pyleoclim/utils/correlation.py
Type:      function

Example case I#

First let’s generate a list of \(p\)-values with size 10 to emulate 10 tests.

np.random.seed(2333)
pval_list = np.random.random(10)/10
print('p-value list:\n', pval_list)

fdr_res = fdr(pval_list, method='original')
print(fdr_res)

p-value list:
 [0.05297776 0.08944875 0.02423382 0.09269607 0.01918608 0.03875847
 0.02873022 0.07296763 0.04393807 0.03032502]
None

In the case above, we get None, which means that there are no significant tests. Note that the 3rd test has p-value equals to 0.024, which will be treated as significant in traditional false positive test where the significance level \(\alpha=0.05\). Now let’s test another case below.

Example case II#

np.random.seed(2333)
pval_list = np.random.random(10)/15
print('p-value list:\n', pval_list, '\n')

fdr_res = fdr(pval_list, method='original')
print(fdr_res, '\n')

print('Index\tp-value')
for idx in fdr_res:
    print(f'{idx}\t{pval_list[idx]}')

p-value list:
 [0.03531851 0.0596325  0.01615588 0.06179738 0.01279072 0.02583898
 0.01915348 0.04864508 0.02929205 0.02021668] 

[2 4 5 6 8 9] 

Index	p-value
2	0.01615588065179053
4	0.01279071873650784
5	0.02583898106901107
6	0.019153480900810392
8	0.029292048677502642
9	0.020216678222839013

In this 2nd case, the fdr function returns a list of indices, which indicates the significant ones among the 10 tests. Note that the 1st test has p-value equals to 0.035, which will be treated as significant in traditional false positive test where the significance level \(\alpha=0.05\).

With the EnsembleSeries.correlation() method, this FDR approach will be applied by default when determine the significance of each p-value.

Correlations and their significance with Pyleoclim

Contents

Correlations and their significance with `Pyleoclim`#

Authors#

Preamble#

Introduction#

Goals#

Keywords#

Pre-requisites#

Relevant Packages#

Data Description and Processing#

Demonstration#

The Corr object#

Pairwise Correlations Three Ways#

Field correlations#

Correlations with `EnsembleSeries` objects#

`EnsembleSeries` vs `EnsembleSeries`#

Takeaways#

Acknowledgments#

About `apply_fdr` in `EnsembleSeries.correlation()`#

Basic idea of FDR#

The underlying FDR function in `Pyleoclim`#

Example case I#

Example case II#

Correlations and their significance with Pyleoclim

Contents

Correlations and their significance with Pyleoclim#

Authors#

Preamble#

Introduction#

Goals#

Keywords#

Pre-requisites#

Relevant Packages#

Data Description and Processing#

Demonstration#

The Corr object#

Pairwise Correlations Three Ways#

Field correlations#

Correlations with EnsembleSeries objects#

EnsembleSeries vs EnsembleSeries#

Takeaways#

Acknowledgments#

About apply_fdr in EnsembleSeries.correlation()#

Basic idea of FDR#

The underlying FDR function in Pyleoclim#

Example case I#

Example case II#

Correlations and their significance with `Pyleoclim`#

Correlations with `EnsembleSeries` objects#

`EnsembleSeries` vs `EnsembleSeries`#

About `apply_fdr` in `EnsembleSeries.correlation()`#

The underlying FDR function in `Pyleoclim`#