Regression spectra (negative values)

[saltwater-tensor]

Hi team,

I was calculating the regression spectra using the function: spectral.regression_spectra()

I am using normalized alpha.

Code

f, psd, coh = spectral.regression_spectra(
    data=data,
    alpha=alpha,
    sampling_frequency=250,
    frequency_range=[1, 100],
    window_length=1000,
    step_size=20,
    n_sub_windows=8,
    return_coef_int=False,
)

Problem

There are negative values in the cpsd.real being passed to the function def coherence_spectra(cpsd, keepdims=False).

https://osl-dynamics.readthedocs.io/en/latest/_modules/osl_dynamics/analysis/spectral.html#coherence_spectra

I don't understand why would we have negative values in the spectrogram that is being calculated using Welch's method.

What is going on here?

[cgohil8]

Normally, we wouldn't use the normalised alphas to calculate the mode-specific spectra. However, I think we standardise the alphas before regressing, so it might not make a difference. Maybe just to double check can you run the function using the raw alphas you get from model.get_alpha()?

[saltwater-tensor]

Yes in the regression spectra code the alphas are normalised.

And yes I can run the function using raw alphas.

[saltwater-tensor]

I ran the scipy.signal spectrogram function on the same dataset. I don't get back any negative values in the spectrogram when I use that function.

[saltwater-tensor]

Just to be clear there are negative values in the actual PSD/auto spectral density which is a problem. The real part of cross spectral density can have negative values which is fine.

[RukuangHuang]

Normally, we wouldn't use the normalised alphas to calculate the mode-specific spectra. However, I think we standardise the alphas before regressing, so it might not make a difference. Maybe just to double check can you run the function using the raw alphas you get from model.get_alpha()?

This might be unrelated but I think the normalisation in normalised alphas is different from the normalisation before regression. The formor ensures they are normalised by the trace of the covariances, but the latter z-transforms the alphas.

[cgohil8]

Sorry for the delay. The model to calculate the mode PSDs is a GLM (General Linear Model):

P_t = alpha_jt P_j + P_0 + eps

where P_t is the spectrogram calculated using Welch's method. alpha_jt are the (non-normalised, i.e. raw) mixing coefficients inferred by DyNeMo (these are obtained using the model.get_alpha() method), P_j are the mode PSD - this is what we're primarily interested in, P_0 is the static (time-averaged) PSD, and eps is the residual. Note, the PSDs are functions of frequency: P_t = P_t(f), P_j = P_j(f), P_0 = P_0(f).

P_t is always positive. This means P_0, which is the time-average of P_t, is also always positive. However, P_j can be negative for certain frequencies which represents for a particular mode then that mode "activates" (increases in mixing coefficient) it results in less power in that particular frequency. In practice however the mode PSDs, P_j, generally are positive for most frequencies and and negative values are very small.

I would advise just plotting P_j (mode PSDs, referred to as the coefficients in the code) and P_j + P_0 (P_0 is referred to as the intercept in the code), do these look reasonable?

Note, if you're passing return_coefs_int=False, then spectral.regression_spectra() should return P_j + P_0 and this should not be negative, maybe there is a bug if this is the case.

[saltwater-tensor]

P_t = Spectrogram for a window designated by t?
P_j = Is this the regression coefficient that scales alpha linearly to produce P_t?

Then when I use return_coefs_int = False will P_j s be omitted? Or just the intercept will be not be returned.

If that is the case then for clarity does the "int" part of return_coefs_int refer to intercept?

So this means that if for example when we use coherence returned from the function regression spectra then it possible that the values are negative and/or greater than 1?

[saltwater-tensor]

Since you mention that P_j can be negative therefore cpsd values can be negative.
And therefore the values that go into the following equation can be negative.

coh[i, j, k] = abs(cpsd[i, j, k]) / np.sqrt(
cpsd[i, j, j].real * cpsd[i, k, k].real

That is why I get warnings and errors for np.sqrt.

Remedy

But this has to be fixed. Can't we simply take the absolute value in order to normalize? Sign has no meaning here for normalization. Maybe directly of the complex number.

Absolute value is being used in the numerator anyway. Although it is being used to calculate the magnitude of complex values.

Some questions

This also brings into question the concept of using regression spectra to calculate coherence or PSD.

What are we assuming here? That it doesn't matter what the sign of regression coefficient is for a specific frequency bin it just depends upon the magnitude?

Example

2 Mode DyNemo model:
For mode 1: 15Hz frequency bin PSD has a negative coefficient (P_j).
Mode 2: 15 Hz frequency bin has a positive P_j (it should right by design?)
But then if we take the absolute value for PSD or coherence in both the cases what does it mean?
Both the modes capture power in the 15Hz frequency?

Then how do we spectrally segregate the modes? Or what is unique about them?
I mean clearly from the regression analysis the sign is proving to be important but then it is does not show
up in the spectral analysis.

Any other interpretation would also be helpful.

[saltwater-tensor]

Sample output

NAN VALUES

:1: RuntimeWarning: invalid value encountered in sqrt

abs(cpsd[i, j, k]) / np.sqrt(cpsd[i, j, j].real * cpsd[i, k, k].real)
:1: RuntimeWarning: invalid value encountered in sqrt
array([9.9997532e-01, 9.9570370e-01, 1.0379384e+00, 1.2386230e+00,
1.3821068e+00, 5.0255984e-01, 2.6649192e-01, nan,
8.8342857e+00, 2.5248156e+00, 1.0680170e+00, 5.0258040e-01,
2.2400749e+00, 1.1725411e+00, 7.1887213e-01, 7.1562386e-01,
1.0495523e+00, nan, 8.7028301e-01, 1.0597816e+00,
2.0021662e-01, nan, 1.1709763e+00, 4.3582582e+00,
1.4349971e+00, 1.0539809e+00, nan, 1.2518290e+00,
3.4332504e+00, 1.3921651e+00, nan, 1.5290527e+00,
1.7810689e+00, nan, nan, 1.0318611e+00,
1.5277568e+00, 2.0876372e+00, 4.9194255e+00, nan,
1.0363775e+00, 4.0008073e+00, 1.0912158e+00, 4.4088712e+00,
3.0298977e+00, nan, nan, 1.3054664e+00,
9.7680277e-01, nan, nan, nan,
6.7534864e-01, 1.3301733e+00, nan, 2.7986636e+00,...

[cgohil8]

Did you bandpass filter the data you pass to spectral.regression_spectra? If so, what range?

[saltwater-tensor]

1-125 Hz

[JieLiang-ustc]

I meet the same problem. I get negative value in psd compute when usd spectral.regression_spectra after training DyNeMo. However, when I train HMM using the same data, I didn't get any negative value.

[cgohil8]

Apologies for the delay on this. Will get back to very soon.