gainFinder cannot regress predicted pupillary response

[nicobast]

The code is now modified to fit our oddball task. However, with sample data the model fails to predict the pupillary response at all. The prediction remains a flat line and likely refers to this warning:

`Warning: X is rank deficient to within machine precision.

In regress (line 84) In gainFinder (line 114) In fitModel (line 38)`

I fear that there might not be enough saccades in our short trials (<0.7s), but any suggestions would be highly appreciated.

[csb0]

That error suggests the regression itself is failing — i.e., there are too many near-zero values in either the generator functions, filter, or trial-averaged pupil responses (data). So I would do a few diagnostic checks. Plot the estimated filter (which the code does automatically I think) and make sure it's not flat or unexpected in some other way. Plot the generator function for each trial type. If either the generator functions or filter have issues you could a get a rank-deficient design matrix in the regression, leading to this error. You should also plot the trial-averaged pupil responses for each trial type that you are fitting to make sure they look as expected.

I can think of other things to try next, if you don't discover the issue by doing the above.

Short trials are not an issue, per se, because you are estimating the saccade rate over many trials.

[nicobast]

Thanks, for your ongoing support. You are right. The estimated filter ist flat (attached). I am not sure, why this is the case, as there is substantial pupillary response and there are relatively frequent saccades. Do you have any suggestions, how to proceed with trouble-shooting?

[nicobast]

and as an addition, also the generator function. The downward outlier do not seem to be very generator-function like?

[csb0]

The trial-averaged pupil response and generator function look fine. Those downward parts are from missing data and don't significantly affect the model fit.

In the figure with the filter that has legends "data" and "parametric fit," the "data" curve (which represents the saccade-locked pupil response) should not be flat. That is likely the main issue. Can you first plot the data curve by itself? In the main script, fitModel.m, it is called d.sacIrf. If that's flat, then there was a problem estimating the saccade-locked pupil response from your data (with deconvolution). Can you check if that's true first?

[nicobast]

Thanks, your are right. The problem seems to be with estimating the linear filter. The IRF from the data (d.sacIrf) seems fine

the problem seems to be occur with this code line

[parametricLinearFilter, params_PuRF, Rsq_PuRF, normFactor_PuRF] = fitParametricPuRFfromData(IrfAvg,d.sampleRate./d.downsampleRate,d.sampleRate,1); which returns the flat line (zeros).

[csb0]

Great!

Because your trials are so short, the saccade-locked pupil response is very distorted (notice the oscillations) and can't be used to parametrically fit the filter. This is actually an expected situation given our theory (and I should add something to the documentation about this and an error-check/message for this case).

The quickest solution is to use a 'default' parametric linear filter, ideally collected in the absence of a task (when the observer just fixates) or with longer trials. I can share ours estimated from 10 observers. It will probably do a decent job of fitting, as there is not a huge amount of inter-observer variability in the filter (we have a SI figure showing the variability among our observers).

If you tell me the sampling rate of your pupil data, I will make this and share with you.

[nicobast]

thanks, for your diagnostic insight. I might suggest another workaround before approaching your path. I would like to try this, to stay as close to your original analysis. So, I looked into the origin of the flat linear filter line. I set the IRFdur back to its default of 4 seconds (which might not make much sense given our trial durations of 600-700ms).

In addition, the flat linear filter seems to be due to the normalizing of the PuRF amplitude. When I comment out the lines in fitParametricPuRFfromData, it will now produce a linear filter and do (currently poor and unoptimized) fit:

%% Normalize the PuRF amplitude just for plotting. we will later renormalize it when we track gain over time % normFactor = abs(min(parametricLinearFilter))./abs(min(irfData)); % parametricLinearFilter = parametricLinearFilter./normFactor; % nov 11, 2020: checked and it barely affects best-fit sigma value (neglibly) to normalize this or not % parametricLinearFilterOrginalSamplingRate = parametricLinearFilterOrginalSamplingRate./normFactor; You already stated in the paper that the normalization had very little effect. Would it be fine to just drop the normalization?

This is the fit output:

[nicobast]

I would also be interested in implementing the default linear filter. Our sampling rate is 300Hz. (We actually have sampling rates of 120Hz and 300Hz across participants, but the sample data I used is 300Hz).

[nicobast]

I think, I am getting there. With some optimized settings, I will get a rather okay fit in the sample data:

But only possible if the normalization is dropped. Do you think that is okay?

[csb0]

I see. The filter should always have a duration of around 4 seconds. It represents the response of the pupillary muscles to a single impulse of internal neural activity. Therefore, this duration depends primarily on the biomechanical properties of the pupil muscles, not on task or behavior-dependent factors. This is also why it's probably OK to use a default one if you need to. The filter you showed now looks fine also.

The normalization matters if you are trying to compare gain (arousal) estimates across observers or across recording sessions within an observer (only applicable when the eye tracker moves wrt the observer, e.g., if they need to adjust their position and re-calibrate). We can probably figure out why it's not working. Can you plot irfData and parametricLinearFilter at line 46 in fitParametricPuRFfromData.m, in two separate plots? The normalization factor is computed as the ratio of the absolute values of the minima of these two functions.

[nicobast]

Yes, it would matter as I would like to compare gain estimates between groups.

Here is the plot of parametricLinearFilter at line 46 in fitParametricPuRFfromData.m:

However, I cannot plot irfData as it is a constant:

see line from fitModel.m: IrfAvg = nanmean(d.sacIrf); % avg. saccade-locked IRF across all runs This IrfAvg is used as input irfData in fitParametricPuRFfromData.m. In addition in fitParametricPuRFfromData.m, this value is subtracted from itself resulting in zero: irfData = irfData-irfData(1);

Any idea, what I get wrong?

[nicobast]

this is the actual d.sacIrf. Did you mean that by "IrfData"?

[csb0]

Ah, that is the problem. irfData is supposed to be a matrix with multiple runs in it. Because you're using just one run, it is averaging over time rather than runs. I will fix this in the code to prevent this from occurring in such a case.

I did so, and the fix was just to change line 32 in fitModel.m to IrfAvg = nanmean(d.sacIrf,1);

[nicobast]

Great, thanks. This seem to have fixed it. I am now able to model it without errors. The fit is now worse, but will try to optimize the parameters for saccade fixation and then come back.

Here is the output:

I will also try to split the 1400 trials into distinct runs and differentiate between trial types.

[csb0]

The saccade-locked pupil response depicted in the figure above (bottom right, black curve) is highly distorted, likely because your trials are quite short. The filter resulting from that probably shouldn't be used to fit the data, and may be why the model fit is so poor. The filter's minimum is abnormally early.

I have provided a filter fit from many observers' data that will likely do a better job. It's now under /old/parametricLinearFilter.mat in the github repo, sampled at 500 Hz. You can use Matlab's built-in resample function so that it is at 300 Hz instead or whatever sampling rate you need. And then just plug that in on line 49 of fitParametricPuRFfromData, in the bold text, so that you normalize it's height for each run or observer separately.

parametricLinearFilter = parametricLinearFilter./normFactor;

Also, is this data from one observer or many? The model works best if used on individual observer's data, one at a time.

[csb0]

There is another alternative, but it involves collecting around 10 min of data from each observer.

[nicobast]

Also, is this data from one observer or many? The model works best if used on individual observer's data, one at a time.

This is from a single sample participant. However, will try to apply the model to a sample of n = 300. I will try to implement your proposed solution and come back to you

[nicobast]

There is another alternative, but it involves collecting around 10 min of data from each observer.

Yeah, this will not be an option as it is post-hoc and the data already collected several years ago.

[nicobast]

% TESTING - used fixed linear filter from Burlingham et al. 2022
S = load('parametricLinearFilter.mat');
parametricLinearFilter_fixed = S.parametricLinearFilter;
parametricLinearFilter_fixed_downsampled = resample(parametricLinearFilter_fixed,3,5);

%% Normalize the PuRF amplitude just for plotting. we will later renormalize it when we track gain over time
normFactor = abs(min(parametricLinearFilter))./abs(min(irfData));
%parametricLinearFilter = parametricLinearFilter./normFactor; % nov 11, 2020: checked and it barely affects best-fit sigma value (neglibly) to normalize this or not
parametricLinearFilterOrginalSamplingRate = parametricLinearFilterOrginalSamplingRate./normFactor;
parametricLinearFilter = parametricLinearFilter_fixed_downsampled./normFactor;`

Based on your suggestion, I made the adaptations above will leads to an Rsq = 0.83 and looks fine. Thanks!

[csb0]

Good to hear!