I have a GLM model with 4 event regressors. 3/4 events occur 7-12 times per scanning session, and one event occurs 1-4 times per scanning session. I have a total of 31 participants.

However, for some participants, the sparse event occurs more frequent than for others, resulting in an unbalanced count both across participants and across conditions.

You can technically analyze the data just as you would with any other GLM, but you just have to keep in mind that the error surrounding the beta coefficients for your sparser event will probably be larger than with your more frequent events. This means that your test statistics (T/Z) for that sparse condition will be likely also be smaller, since these variables take variance of beta coefficients into account. However, you usually aren’t converting to T/Z stats until the end (e.g., after contrasting betas from different conditions).

Hey thanks for the rapid feedback,
I performed my analysis, which yielded some very interesting results, but I was given some feedback that because there is the problem of heteroskedasticity, the assumptions of the traditional summary statistics approach are not valid. Therefore, I cannot move the betas at the second level and do inference on them using a t-test for example.
Because of that, I am trying to find additional ways to do inference. Is there another approach, or a way to justify the summary statistics approach ?

I think that in your case, a paired t-test is valid. It may be suboptimal in terms of power, but I don’t see a fundamental for it to be invalid.
By paired t-test, I means that you want to test the reponse difference btw the target condition with fewer trials and another one.
HTH,
Bertrand Thirion