FOOOF aperiodic (over)estimation

[stefandukich]

Hello,

I am working with the FOOOF method using the Brainstorm implementation (FOOOF_matlab) on resting-state data (6 minutes total, 2 s epochs, 50% overlap, 0.3 Hz highpass filter).

I have noticed that the method often inaccurately estimates the spectral peaks (e.g., identifying one large peak instead of two, or completely missing a peak). To address this, I thought of estimating the periodic power manually by calculating the difference between the raw power and the estimated aperiodic component. When I do this, however, I get negative values at the "edge" frequencies (see the attached figure; averages across 70 participants).

I have tried using the default parameters as well as varying them (e.g., narrowing down the freq. span), but the overestimation of the aperiodic component persists. Also, I am using the 'fixed' aperiodic mode as this is what is typically used (which might be a questionable choice, of course).

I would be very keen to hear your opinion on this. Papers often do not plot the raw spectra versus FOOOF estimates, so I am unsure if this is a common occurrence or something specific to my data. Any advice would be appreciated.

Thanks, Stefan

[Raymundo.Cassani]

@stefandukich

I have noticed that the method often inaccurately estimates the spectral peaks (e.g., identifying one large peak instead of two, or completely missing a peak).

Have you taken a look to the model-selection addition to the spectral parametrization method?
It aims to improve the peak detection and reduce the subjective parameter choices.

https://neuroimage.usc.edu/brainstorm/Tutorials/Fooof#Model-selection_.28ms-specparam.29

@stefandukich

When I do this, however, I get negative values at the "edge" frequencies (see the attached figure; averages across 70 participants).

This is quite common, and rather than on the edge frequencies, it appears on the tails of the (Gaussian) peaks, thus it is also common to find it in the middle of the PSD in between peaks, as in this example, the green line is the raw PSD.

---

@Luc, do you want to add something else to the conversation?

[stefandukich]

Hi,

Appogies, it took me a bit to get back to this analysis properly.

First of all, thank you very much for the response!

Model-selection method

I had noticed this method, but I had not looked into it until now. I have now applied it to my data and found that it did improve the estimates. Nevertheless, negative values are still present at lower frequencies (which may be a result of the knee effect) and, depending on the cohort, they appear at higher frequencies as well (potentially due to flooring caused by EMG). Here are the new estimates from the two groups:

Peak estimation issue

When examining individual estimates, I have noticed that the method relatively often identifies a single large peak even when several are clearly present. As my setting is peak_width_limits = [1, 10], this is very likely caused by the higher upper limit (i.e., 10 Hz).

To mitigate this, I could consider lowering the limit to, for instance, 8 Hz… but, I believe that it is good to have the option to fit broader peaks as well. So, I thought to add the following logic to my code (after line 276):

    % Unusual to have just one peak
    % If only one peak is detected, it is typically not correct
    if size(est_pars, 1) == 1
        % disp(chan);
        delta_lim = 0;
        % fprintf('start: %d\n', size(est_pars, 1));
        while true
            delta_lim = delta_lim + 1;
            % peak_threshold = 1.5;
            peak_threshold = opt.peak_threshold;

            % Calculate the effective upper limit for this iteration
            current_upper_lim = opt.peak_width_limits(2) - delta_lim;

            [est_pars, peak_function] = est_peaks(fs, flat_spec, opt.max_peaks, peak_threshold, opt.min_peak_height, ...
                [opt.peak_width_limits(1) current_upper_lim]/2, opt.proximity_threshold, opt.border_threshold, opt.peak_type);

            % Stop condition:
            % 1. If we have found more than 1 peak OR
            % 2. If the upper width limit has dropped too low
            % fprintf('-> %d\n', size(est_pars, 1));
            if size(est_pars, 1) > 1 || current_upper_lim == 4
                break;
            end
        end
    end

In essence, this snippet gives the method an opportunity to fit more peaks by iteratively lowering the upper limit. While I understand that there may genuinely be only one peak in data, I also wonder how often that is actually the case in practice. After all, the model-selection apporach can determine the best number of peaks in the subsequent steps of the original code.

Anyway… This idea helped, but it did not always resolve the issue completely:

Non-linearity in the log-log space

Finally, I have noticed another problem where certain electrodes are simply not linear in the specified frequency range. I am not entirely sure how to deal with these situations, especially as they occur occasionally and inconsistently across different participants and electrodes (BTW, this does not always happen in the peripheral electrodes where you would expect EMG/EOG activity to cause non-linearity).

While I could restrict the range to, for instance, 3-30 Hz, I am unsure if this is the best thing to do. First, this would sacrifice frequency ranges that carry important information. Second, because these frequencies are within the range where true oscillations can be present, the method would not fit the correct model if the limits happen to be where the peaks are in the data.

So, I wonder if there is a good solution to handle these occasional incorrect estimates, or if it is better to treat these data points simply as missing data?

Thank you very much!

Best,
Stefan