GLMSingle failing at point of running type C model (GLMDenoise) with "Matrix is close to singular or badly scaled"

[cracey]

I'm running glmsingle on data I have previously analysed with glmdenoise.
Model typeA (ONOFF) and type B (FITHRF) seem run successfully but while running type B I get the following warning on each 'chunk'...

working on chunk 9 of 9.
Warning: One or more regressors are all zeros; we will estimate a 0 weight for those regressors.
Warning: One or more regressors are all zeros; we will estimate a 0 weight for those regressors.
Warning: One or more regressors are all zeros; we will estimate a 0 weight for those regressors.

Then when switching to running model C it fails at the initial stages with the following error...

*** Saving results to /home/chris/cisc1/projects/forster_isp/firstlevel/CISC12263/GLMsingle/TYPEB_FITHRF.mat. ***
*** DETERMINING GLMDENOISE REGRESSORS ***
*** CROSS-VALIDATING DIFFERENT NUMBERS OF REGRESSORS ***
working on chunk 1 of 9.
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  4.825304e-08. 
Error using olsmatrix2 (line 68)
Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  4.825304e-08.

Error in GLMestimatemodel>fitmodel_helper (line 802)
    f = mtimescell(olsmatrix2(cat(1,design{:}),opt.lambda),data2);  % conditions x voxels

Error in GLMestimatemodel (line 412)
    fitmodel_helper(design,data2,tr,hrfmodel,hrfknobs, ...

Error in GLMestimatesingletrial (line 760)
        [results0(ll+1),cache] = GLMestimatemodel(designSINGLE,data0, ...

Error in glm_singletest (line 106)
    [results] = GLMestimatesingletrial(design,data,stimdur,tr,[outputdir '/GLMsingle'],opt);

The problem appears to be to do with the internal single trial design matrix with all runs combined.

The design matrix (X) that is entered into the ridge regression on line 67 of olsmatrix2.m isn't quite rank deficient (i.e., rank(X)==size(X,2)), yet it is close to singular  - Please see the attached histogram of singular values (returned by an SVD on X) which shows the smallest singular value to be 0.001735.

My design matrix and data both work fine with GLMDenoise, and the results outputs of the FITHRF model look sensible (in both cases, no output beta == 0).
Is there a way to resolve this or is our design just not compatible with glmsingle? Any advice would be greatly appreciated.

My design matrix is available here... 

cr_designmat.zip

and my data are shaped like this... (with a stimdur=1 and tr=0.8)
data =
  1×5 cell array
    {407722×475 single}    {407722×520 single}    {407722×519 single}    {407722×520 single}    {407722×520 single}

Best wishes,
Chris

[kendrickkay]

Hi Chris,

I took a look at the design matrix and was able to reproduce the crash. Which is a good thing.

The fundamental issue here lies, as you suspect, with the design matrix. Here's a quick look:

[kendrickkay]

And if we simply sum across the columns, we see this:

[kendrickkay]

So what this is showing is that you have indicated the fact that a new trial occurs, occasionally, at every successive TR (within a given time window). For example, TRs 6-10 in the attached figure. This is a bit "insane", in the sense that you are trying to estimate a separate BOLD response amplitude every 0.8 s (since your TR is 0.8 s). This causes extremely highly correlated single trial regressors, like:

[kendrickkay]

Is this your intention? Is there an alternative way of coding your design matrix? For example, I noticed in columns 13 14 and 15 there is a string of ones. Did you intend this? Could we substitute a single "onset" for the string of 5 ones? Or maybe something like [T1 0 T2 0 T3] or something like that?

Also, is the design matrix in general matching your intent? Is there really a new 1-s trial that occurs every 2.4 s? (Like 3 TRs of 0.8-s)? So you have 1-s trial and then 1.4-s gap and then the next trial?

If that is correct, that is quite aggressive and fast. GLMsingle should be able to cope, but things get more challenging to estimate once you push things that quickly.

Finally, the larger picture issue, which, arguably is the source of the crash that you are observing is that it appears you essentially have no "down time" during your experiment. Is that correct? If you have close to no dead time, this is going to cause collinearity issues in the sense that we cannot reliably estimate baseline (e.g. through the polynomial regressors). Or maybe can we treat whatever is happening in columns 13, 14, and 15 as reflecting baseline time, so that what you are trying to estimate in the GLM is the BOLD response amplitude relative to that time?

[cracey]

Hi Kendrick,

Many thanks for your detailed response. You are correct there are very few blank periods in this experiment. I had forgotten about this limitation as it's a study scanned some time ago. This was an oversight we made at the design phase but I thought we would be able to work around it.

The trials in columns 13-15 are actually regressors of no-interest. They are the instruction events which runs for 5 seconds at the start of each block of trials. I coded them as consecutive 1s trials to account for the need for a single stimdur with glmsingle. I included them in the model because they are a visual stimulus with some differences from fixation, I think you are right that might be best to leave these unmodeled and include them in the baseline estimation.

I have tried removing these columns from the design and we are left with estimating 648 trials instead of 1188. It doesn't resolve the basic issue (although it does improve things slightly).

Error using olsmatrix2 (line 68)
Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  9.480937e-08.
(previously RCOND =  4.825304e-08.)

I'm not sure if there are any other modifications I can make to the design matrix. The remaining trials are visual or auditory stimulus presentations. Previously when analysing this study using glmdenoise I used a condition split approach to model 4 instances of each trial type. This wasn't really enough for our planned analyses, so I hoped GLMsingle and it's various benefits would improve matters, but full single trial estimates (648) is probably more than we really need also. Is there a way with glmsingle to estimate some amount less than the full 648? And would doing so mitigate this collinearity issue?

Best wishes,
Chris

[kendrickkay]

Aha, I found the problem. Your design goes to the very end of the run! so for the very last nominal trial, you actually, due to the HRF delay, have no chance of estimating the response. What you need to do is to "delete" the coding of the very last few trials (and compensate for it accordingly in your post hoc analyses).

[kendrickkay]

Notice in the above figure the column that looks all blue. I'm pretty sure this is your "last trial" of the first run. This is bad news.

The tool should be changed to warn users about this; I think we've seen this pernicious issue in other designs before. Thanks for flagging it.

[cracey]

Hi Kendrick,

Thanks for figuring this out. I was worried our design was going to turn out to be incompatible with GLMsingle. It now works and is producing sensible results, I'm excited to explore the benefits of the toolbox.

Best,
Chris

[kendrickkay]

Great - hope things seem to process well for your dataset. Let me know if weird issues arise.

[hJosiane]

Hello

I am running into this exact same issue. I get the TYPEA and TYPEB output. when I get to glmdenoise it doesnt work. I deleted the last trials on my design matrix, but the issue is still unresolved. I have read this issue through and through but still haven't found a way to solve it. I might be missing something or have a design that is not GLMsingle compatible.

This is my design matrix
Vdesign.mat.zip

It has 4 conditions which are visual stimuli presented on different frequencies

Please advise

Best,
Josiane

[kendrickkay]

Josiane, there is a "1" in the last row of the first column. Did you remove that "1"? What error message are you seeing? And is your experiment consisting of only 1 run? And are you running the matlab or python version?

[hJosiane]

No, I did not remove that 1. These are the onset times of each stimulus. They have a different duration from 3-6 s so I set stimdur = 4.5. And I am implementing this in Matlab sorry I forgot to mention that. I successfully used GLMsingle on the auditory version of this experiment. These are visual stimuli.

% v1 onset
event1_onset = [1; 34; 59; 85; 112; 143; 174; 236; 253; 259; 298; 307];

% v2onset
event2_onset = [43; 53; 80; 92; 107; 122; 140; 184; 187; 196; 231; 289];

% v3onset
event3_onset = [4; 68; 96; 101; 154; 159; 192; 210; 218; 246; 250; 264; 268];

% v4onset
event4_onset = [9; 13; 64; 130; 136; 149; 168; 202; 207; 214; 222; 282; 286]

[hJosiane]

This is the exact error I am getting:

Error using assert
Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.223027e-09.

Error in olsmatrix2 (line 68)
assert(isempty(lastwarn),lastwarn);

Error in GLMestimatemodel>fitmodel_helper (line 820)
f = mtimescell(olsmatrix2(cat(1,design{:}),opt.lambda),data2); % conditions x voxels

Error in GLMestimatemodel (line 417)
fitmodel_helper(design,data2,tr,hrfmodel,hrfknobs, ...

Error in GLMestimatesingletrial (line 830)
[results0(ll+1),cache] = GLMestimatemodel(designSINGLE,data0, ...

[kendrickkay]

I see. Can you remove that "1" (by setting it to 0)?

(I presume that when you ran the code you got a little warning about "You have specified trial onsets that occur less than 8 seconds from the end of at least one of the runs"...)

[hJosiane]

I initially had it set to 0 and I changed it to 1 to see what happens. I still get the error the same error. The warning I am getting is : Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.223027e-09.

In olsmatrix2 (line 67)

followed by the error
In GLMestimatemodel>fitmodel_helper (line 820)
In GLMestimatemodel (line 417)
In GLMestimatesingletrial (line 830)
Error using assert
Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.223027e-09.

Error in olsmatrix2 (line 68)
assert(isempty(lastwarn),lastwarn);

Error in GLMestimatemodel>fitmodel_helper (line 820)
f = mtimescell(olsmatrix2(cat(1,design{:}),opt.lambda),data2); % conditions x voxels

Error in GLMestimatemodel (line 417)
fitmodel_helper(design,data2,tr,hrfmodel,hrfknobs, ...

Error in GLMestimatesingletrial (line 830)
[results0(ll+1),cache] = GLMestimatemodel(designSINGLE,data0, ...

[kendrickkay]

OK. Could you please send me the full set of design matrices that you are using (for all the runs)?

[hJosiane]

I might be missing your question but I am guessing you asked the one I am using the visual and auditory. The one I sent earlier is for visual (Vdesign); the second matrix I am attaching is the auditory one (Adesign).

Vdesign.mat.zip

Adesign.mat.zip

[kendrickkay]

I noticed that in Adesign, you seem to have a trial that is simultaneously coded as belonging to two different conditions (around row 1 and row 85). This might be the problem. Is it possible for you to remove one of the pair in each of the two cases? The reason that this is a problem is that eventually GLMsingle attempts to create a single trial design matrix where every column corresponds to a distinct trial. I'm not sure how GLMsingle is processing your case but it might be that it ends up with two columns that are literally identical, leading to the singular matrix error

[hJosiane]

yes, it is possible to remove a pair and model fewer trials. I am wondering though would be the case if I am currently running GLMsingle with Vdesign separately. I already modeled trials using the Adesign. Now I am getting the error when I use Vdesign.

…

On Thu, 30 Nov 2023 at 21:02, Kendrick Kay ***@***.***> wrote: I noticed that in Adesign, you seem to have a trial that is simultaneously coded as belonging to two different conditions (around row 1 and row 85). This might be the problem. Is it possible for you to remove one of the pair in each of the two cases? The reason that this is a problem is that eventually GLMsingle attempts to create a single trial design matrix where every column corresponds to a distinct trial. I'm not sure how GLMsingle is processing your case but it might be that it ends up with two columns that are literally identical, leading to the singular matrix error image.png (view on web) <https://github.com/cvnlab/GLMsingle/assets/2728338/02e009fb-a49f-4132-bd6b-dfa9fde4940f> — Reply to this email directly, view it on GitHub <#52 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BELVYSLNCBQ4Y6CSMZEPPETYHE3FFAVCNFSM5UQNOKK2U5DIOJSWCZC7NNSXTOKENFZWG5LTONUW63SDN5WW2ZLOOQ5TONZSGQYTCOA> . You are receiving this because you commented.Message ID: ***@***.***>

[kendrickkay]

OK, I am confused. To help troubleshoot, please prepare a sample call to GLMestimatesingletrial.m that is giving you an error. I probably don't need the variable, but I will need the and and variables. Can you supply this?

[hJosiane]

I am so sorry for the confusion. The experiment had two parts visual stimuli and auditory stimuli. I successfully used GLMsingle for the auditory stimuli. I used GLMsingle with the Adesign, I attached earlier, It ran with no problem. Now the part I am having the error with is the visual part. Here I am trying to model 4 flickering checker-board presented at different frequencies. The design matrix for this is what I am attaching; Vdesign.mat. So these two parts of the experiment are being analyzed independently. You mentioned you will need the.. but you left it out. Did you mean the dataset? I am not sure what you meant so I will provide more context for all the variables I am using to call GLMestimatesingletrial.m I am attaching the input dataset filtered_func_data_S_2mm.nii.gz (after FSL preprocessing , the end goal is to run MVPA after GLMsingle) stimdur = 4.5 opt = struct('wantmemoryoutputs',[1 1 1 1]); [results,designSINGLE] = GLMestimatesingletrial(design.sparseMatrix,img, stimdur,1.5,[outpath '/GLMsingle'],opt); I am attaching both the design, the data,the exact script, and the exact error I am getting. You don't need the whole script, the first part is only looping through all the runs and subjects I have but I attached too because you can see how I am calling GLMestimatesingletrial.m more clearly from there. filtered_func_data_S_2mm.nii.gz <https://drive.google.com/file/d/15fPBJZMi6XwrK2ghzcMK0GCg-P1xYoHV/view?usp=drive_web> I truly appreciate you for taking your time to help me with this issue. I can't thank you enough. I am happy to provide any more details if needed. Thank you so much

…

On Thu, 30 Nov 2023 at 22:19, Kendrick Kay ***@***.***> wrote: OK, I am confused. To help troubleshoot, please prepare a sample call to GLMestimatesingletrial.m that is giving you an error. I probably don't need the variable, but I will need the and and variables. Can you supply this? — Reply to this email directly, view it on GitHub <#52 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BELVYSO66EBAX7J6RMB33XTYHFEFTAVCNFSM5UQNOKK2U5DIOJSWCZC7NNSXTOKENFZWG5LTONUW63SDN5WW2ZLOOQ5TONZSGQ2TQOI> . You are receiving this because you commented.Message ID: ***@***.***>

[kendrickkay]

Hello. Sorry for the long delay, was dealing with a lot of fires.

I took a look. It appears I only got the .nii.gz file in your message. (Can you send the other things, script, design, etc.)?

In considering your message, a few thoughts:

• You can consider simultaneously modeling both the auditory and visual parts of your experiment in a single call to GLMestimatesingletrial.m. Like, give all the data to that function and prepare appropriate design matrices. There is some reason that doing it simultaneously (for the data from a scan session) will actually give better results.
• The .nii.gz file you sent has 208 volumes at a TR of 1.5 s. The data look fine. The data look somewhat smooth (perhaps you spatially smoothed a little?), and it looks like you applied a brain mask (this is fine from GLMsingle's point of view).
• However, the Vdesign.mat file you provided earlier has 298 rows (it is 298 x 4). The number of rows in the design matrix should match the number of fMRI volumes you prepare. Are the data and design reflecting the same TR (e.g. 1.5 s)?

If you can supply the full context and error message, I can take a look (more promptly this time).

[hJosiane]

Hello No problem I also apologize for not replying to your email earlier. I have been caught up in deadlines. I do still want to use GLMsingle and still have the same error so thank you for the email. Those comments are very useful: 1. I will try to make a design that combines both matrices to get better results 2. yes, prior to GLMsingle we did smoothing and warped the data in standard space. The end goal is to do multivariate pattern analysis and we wanted to be able to do between-subject analysis, thats why we did this so that we can just work in standard space. 3. Yes, the Vdesign has 298 rows and the volumes are 208 (TR 1.5S). The visual stimuli weren't the main part so they were displayed longer. I can simply remove the last blocks in the design Lastly, I apologize for not attaching the code I mentioned and the design. I am attaching them here I will attach the dataset for one subject again. Script is GLMsingleAll_Visual.m The first design is Vdesign.mat ( V stands for visual stimuli), and the second design is Adesign. (A stands for auditory stimuli). This is the dataset again filtered_func_data_S_2mm.nii.gz.zip <https://drive.google.com/file/d/1qD1YucGSMuyXZ6-ZSbXVAduF5Syx5YN1/view?usp=drive_web> Thanks again

…

On Wed, 6 Dec 2023 at 09:05, Kendrick Kay ***@***.***> wrote: Hello. Sorry for the long delay, was dealing with a lot of fires. I took a look. It appears I only got the .nii.gz file in your message. (Can you send the other things, script, design, etc.)? In considering your message, a few thoughts: - You can consider simultaneously modeling both the auditory and visual parts of your experiment in a single call to GLMestimatesingletrial.m. Like, give all the data to that function and prepare appropriate design matrices. There is some reason that doing it simultaneously (for the data from a scan session) will actually give better results. - The .nii.gz file you sent has 208 volumes at a TR of 1.5 s. The data look fine. The data look somewhat smooth (perhaps you spatially smoothed a little?), and it looks like you applied a brain mask (this is fine from GLMsingle's point of view). - However, the Vdesign.mat file you provided earlier has 298 rows (it is 298 x 4). The number of rows in the design matrix should match the number of fMRI volumes you prepare. Are the data and design reflecting the same TR (e.g. 1.5 s)? If you can supply the full context and error message, I can take a look (more promptly this time). — Reply to this email directly, view it on GitHub <#52 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BELVYSNQP4O6YWD6Q2C4KZTYIB3SRAVCNFSM5UQNOKK2U5DIOJSWCZC7NNSXTOKENFZWG5LTONUW63SDN5WW2ZLOOQ5TONZXGY3TCOA> . You are receiving this because you commented.Message ID: ***@***.***>

[kendrickkay]

Hi - I see a link for the dataset (filtered_func_data_S_2mm.nii.gz.zip). But I do not see the other things you mentioned (e.g. GLMsingleAll_Visual.m; Vdesign.mat).

[hJosiane]

They are attached to the email. There are 4 attachments to the previous email. [image: Screenshot 2023-12-13 at 12.12.37 PM.png] Best regards, Josiane

…

On Wed, 13 Dec 2023 at 10:39, Kendrick Kay ***@***.***> wrote: Hi - I see a link for the dataset (filtered_func_data_S_2mm.nii.gz.zip). But I do not see the other things you mentioned (e.g. GLMsingleAll_Visual.m; Vdesign.mat). — Reply to this email directly, view it on GitHub <#52 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BELVYSKREGGXIYZ5YBKJLE3YJHDZJAVCNFSM5UQNOKK2U5DIOJSWCZC7NNSXTOKENFZWG5LTONUW63SDN5WW2ZLOOQ5TOOBUGM4TQMA> . You are receiving this because you commented.Message ID: ***@***.***>

[kendrickkay]

Sorry - It appears github is stripping the attachments. You can just send directly to ***@***.*** ***@***.***>

…

On Dec 13, 2023, at 11:13 AM, hJosiane ***@***.***> wrote: They are attached to the email. There are 4 attachments to the previous email. [image: Screenshot 2023-12-13 at 12.12.37 PM.png] Best regards, Josiane On Wed, 13 Dec 2023 at 10:39, Kendrick Kay ***@***.***> wrote: > Hi - I see a link for the dataset (filtered_func_data_S_2mm.nii.gz.zip). > But I do not see the other things you mentioned (e.g. > GLMsingleAll_Visual.m; Vdesign.mat). > > — > Reply to this email directly, view it on GitHub > <#52 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/BELVYSKREGGXIYZ5YBKJLE3YJHDZJAVCNFSM5UQNOKK2U5DIOJSWCZC7NNSXTOKENFZWG5LTONUW63SDN5WW2ZLOOQ5TOOBUGM4TQMA> > . > You are receiving this because you commented.Message ID: > ***@***.***> > — Reply to this email directly, view it on GitHub <#52 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAU2DEWUJ4AS4BGWJSCHCSLYJHO3DAVCNFSM5UQNOKK2U5DIOJSWCZC7NNSXTOKENFZWG5LTONUW63SDN5WW2ZLOOQ5TOOBUGUYTKMI>. You are receiving this because you commented.

[kendrickkay]

Sorry - It appears github is stripping the attachments. You can just send directly to ***@***.*** ***@***.***>

Whoops - looks like the email address is masked too. You can just find my email online and email me directly.

[kendrickkay]

Hi Josiane,

I took a look at your materials. I have a few observations.

The data (filtered_func_data_S_2mm.nii) has 208 volumes. With a TR specified at 1.5 s, this means your fMRI data are 312 s long.

The design matrix (design.sparseMatrix) is 298 x 4. The intention of GLMsingle is that the rows of your design matrix should correspond to the TR. Hence, there should be 208 rows in your design matrix (so that it matches your fMRI data).

The reason that it is crashing is that GLMsingle automatically drops the extra rows of the design matrix (i.e. it keeps only the first 208 rows). After doing that, there are trials that are dangerously close to the end of the run, and so that causes the "singular error".

Another observation is that calling GLMsingle with only one run means that you will not really be able to take advantage of the GLMdenoise and RR components. So if you have multiple runs, you should call GLMsingle with all of the data.

See if this makes sense. If you are still getting errors, please open a new github Discussions thread.