T_period in method "seasonal"

[ClarGl]

Hello, thank you for the excellent rEDM package, which works really well for addressing my research question. However, I came across an issue where I was hoping to get your help: What is the temporal resolution of “T_period” when creating a matrix with surrogates and applying the “seasonal” method? I assume that the default setting is a value of one for T_period. However, I don’t understand whether this value corresponds to one week, month, day, or year—especially since the first column of my data is a string and does not include any date.

Additionally, I am studying data where I expect external forcing due to solar radiation (diurnal fluctuations) rather than seasonal external forcing (seasonality). Is there a way to test for this external forcing on a daily basis, particularly if I can mathematically approximate this diurnal, external forcing using the same spline function applied for seasonality?

Thanks for your help!

[SoftwareLiteracy]

Hi,

It looks to me like the SurrogateData(method = "seasonal") function is broken.

For example:

library(rEDM)
N = 300
t = seq(1,N)
ts = sin( 5*pi*t/N ) + rnorm( N, 0, 0.1 )
S = SurrogateData( ts, method = "seasonal", num_surr = 3, T_period = 1, alpha = 0 )
Error in smooth.spline(c(I_season - T_period, I_season, I_season + T_period),  : 
  need at least four unique 'x' values

I believe this function came from the legacy code rEDM 0.7, but apparently not functionally validated. Looking at the SurrogateData() code for method = "seasonal" it makes no sense to me. This looks to be a handcrafted algorithm for a specific use-case rather than a generalized function.

I'm not clear on the application you describe. It sounds like you have data and you want to assess impact of insolation. The sample rate and timescale of your data will be important in determining how to assess insolation. Generating surrogates can then be done based on whether you want to assess diurnal fluctuations or changes in the diurnal pattern itself.

fyi: The python package pyEDM has a SurrogateData() function for seasonal data. It fits a spline to the data, adds resampled (shuffled) residuals to the spline, with optional Gaussian noise. That should be easy to code in R.

[ClarGl]

The SurrogateData function works with my data, so I don't believe it is broken. However, I couldn't find any information about the temporal resolution of T_period, an input parameter for this function.

[SoftwareLiteracy]

If you are happy with the function that is great.

As mentioned, I do not understand the code:

    n = length(ts)
    I_season = suppressWarnings( matrix( 1:T_period, nrow = n, ncol = 1 ) )
    
    # Calculate seasonal cycle using smooth.spline
    seasonal_F = smooth.spline(
      c(I_season - T_period, I_season, I_season + T_period), c(ts, ts, ts) )
    
    seasonal_cyc   = predict( seasonal_F, I_season ) $ y
    seasonal_resid = ts - seasonal_cyc
    
    return(sapply(1:num_surr, function(i) {
      seasonal_cyc + sample(seasonal_resid, n) + rnorm(n, 0, alpha)
    }))

It would be great if you can decipher the intent/meaning of T_period and answer the question here for others.

[nonlinearnature]

This should be the seasonal surrogate method described in https://doi.org/10.1073/pnas.1607747113, which calculates a spline-smoothed average, calculates the deviations, then randomizing the temporal order of those before adding them back to the smooth average.

In this case, T_period is the period of the annual cycle in time-series steps. If it is monthly data, T_period should be 12. If the input time series is not regularly sampled, you will have to be a bit more clever!

[SoftwareLiteracy]

Indeed, the prescriptive algorithm is clear and sensible 👍. The code looks to have some extra sauce with a pre/append of ts.

In any case, the definition/usage of T_period in the docs perhaps deserves tweaking. Help offers:

T_period: the period of seasonality for seasonal surrogates (ignored for other methods)

[ClarGl]

I agree that the documentation for T_period lacks clarity. For instance, it’s unclear whether T_period should be set to 365 when working with daily resolution data?

Thank you for your assistance nonetheless!

[nonlinearnature]

It depends on what your null hypothesis is. If you are seeking to reject the null hypothesis that annual climate forcing is creating prediction / cross-prediction skill in time series that otherwise lack dynamics / dynamic interactions, then yes. However, depending on the type of data, e.g. daily maximum chlorophyll in an estuary, a more important null test might be against the changes in tides over the period of lunar declination (27 days).

[ClarGl]

That makes sense to me—thank you for the example with the different null hypotheses. Could you provide the code for testing daily maximum chlorophyll in relation to changes in the tide? Alternatively, could you recommend a publication where I can find the appropriate code? This would help me set up the code for my null hypothesis, as I'm having difficulty interpreting T_period from the R documentation. Thanks again!