More Generic API for GLM

[Edouard2laire]

Hi @thomas-vincent

This discussion is related to the last lab meeting we had with the NIRSTORM team and #117.

Context:

• We want a process to enhance the linear detrending that offer brainstorm by also removing the slow frequency (eg instead of using a high-pass filler)
• As you explained in #115 , we want a process to remove physiological fluctuation coming from superficial layer
• We want, in the GLM, to be able to add easily new regressor such as the short-separation channel (what you added in #115)

Those 3 points are all related to GLM. In the first two we say Y=XB + e and we want to return e, therefore removing the effect of the regressor in X. In the first case X contains DCT, in the second X contains the value of the SSC. The thirds case, is also similar but we want B.

Problem:

As you have seen, adding a new regressor in X is not easy and therefore your addition in #115 is really hard to maintain and to understand for exemple if we were to add new regressor such as accelerometer data.

What I suggest:

Having a small structure called model that contains our GLM, and a set of function that will interact with it to do what we want : - add regressor, fit the data, display the design matrix...

The structure:

• X : design matrix (time*n_regressor)
• reg_names: 1* n_regressor cell containing the regressors name.
• n_roi : Number of regressor of interest. (present in the first column of X)
• ntime : number of time point
• fs : signal frequency

API:

• model=initialize_model(ntime,fs) return an empty structure

• model = add_regressors(Model, Regressor_type, params) Add regressors to the design matrix. Call one of the following function depending on the regressor type (eg "event", "channel", "constant", "DCT"). Can be extended easily to have more regressor type. model = add_regressors(Model, "event", Events, basis_choice, hrf_duration) model = add_regressors(Model, "channel", channel_names, [convolve_with_hrf=0, basis_choice, hrf_duration]) Note: if convolve_with_hrf=1, then convolute the regressor with HRF like for event based regressor. Typically if the channel is AUX that contain the paradigm of the experiment model = add_regressors(Model, "constant") Add a regressor that just contains 1. model = add_regressors(Model, "DCT", frequences, names)

• display_model(model) Display the model

• [B, proj_X] = model_fit_B(model, y, method='SVD') Return the B, and the corresponding projector base on the method. Method is either 'SVD' or 'pinv' In case of pinv, proj_X = ( X^{T}X)^{-1} X^{T}

• [B, covB, dfe, residuals, mse_residuals] = model_fit(model, y, [method='precoloring', fitoptions]) Call the corresponding method to fit the model depending on the chosen method. ( pre-coloring by default)

• model = model_apply_IIR_filter(model, filter_type, low, high) Apply the IIR filter on the regressor of interest (eg the n_roi first column of X)

Its quite easy to do, as all the material exist. but can really simplify the code for case 1, and 2 : 1/

model=initialize_model(ntime,fs) 
model= add_regressors(model, "DCT",  freq, freq_name);
B=model_fit_B(model, y);
y= y - X*B;   

2/

model=initialize_model(ntime,fs) 
model= add_regressors(model, "channel", superfical_chann);
B=model_fit_B(model, y);
y= y - X*B;   

Don't hesitate if you have comment :) Best regards, Edouard

[Edouard2laire]

Hi @thomas-vincent

• Detrending

This is working nicely, and can remove trend in the data and is able to keep the mean if applied on raw data.

• Remove superficials noise This remove superficial noise; by adding the timecourse of the corresponding channels in the design matrix. Channels are either detected by distance or by manual entry. I only add HbO and HbR data as HbT=HbO+HbR and that might cause problem when estimating B

Maybe we want to add a condition to check that the regressor at not too correlated and then maybe add only meaningful time-course (maybe add the mean, or do more advance stuff like PCA..)

• GLM

The main changes is that before, when applying the high-pass filter on the design matrix. We were assuming that the constant regressor was at the end, or with your changes in #115 , It was no longer true and it might have cause trouble. So I fixed that.

I also add the possibility to let the user say which filter he used for detrending (either IIR, or the new detrending process that I would recommend)

I therefore have a question when filtering the design matrix. Let say we have this simple design :

When using IIR filter on it (low-cutoff = 0.05 = 1/200), we obtains :

When using detrending ( only constant and linear), we obtains :

When using detrending ( only constant and linear and DCT cutoff = 0.05 = 1/200 ), we obtains :

I have the feeling that we are introducing distortion to the HRF when using IIR filter, or when adding DCT in the detrending. Maybe we want to only remove the mean ? But I guess it's mandatory to apply the same filter that was applied on the data on the design matrix. What do you think ?

Another change is that I add the possibility to add nuisance regressor such as short-distance separation channel. We have here the same issue as before as adding more regressor might create trouble of conditioning of the matrix. What do you think of that ? I saw in a recent article about connectivity that they are only regressing out the mean of the short-distance channel. Might be a good idea. What do you think ?

Best regards, Edouard

[Edouard2laire]