FIR analysis with variable intervals

there is variable time interval both before and after Y, which gives me the impression that I will need to build a custom design matrix for each participants each run. Is there anything I should be aware of when running 3dDeconvolve?

To clarify, are you trying to capture the BOLD responses for the four intervals of X1-Y1,Y1-Z1, X2-Y2, and Y2-Z2? How many trials do you have for each of those four intervals? Does the duration vary across trials for each specific interval?

To complicate matters further, not all participants completed the same amount of runs so this variability should be incorporated in the second-level estimates. I believe this means that I will need to run 3dLMER on the second level correct?

If you want to estimate the HRF at the individual level, 3dMSS would be a suitable option as discussed here.

What is the correct procedure to compare my two conditions FIR time courses? Can one do non-parametric cluster correction for each beta of each TR.. and then somehow correct for multiple comparisons additionally due to the multiple TRs?

It would be better to treat each HRF as a whole continuous curve rather than a number of discrete data points. See the same 3dMSS discussion mentioned above.

Gang