Convenient way to manipulate X.xmat.1D design matrix for 3dREMLfit refitting

AFNI Message Board
,
1

Hello,

I have five runs of task fMRI data. They are processed by afni_proc.py. I would like to compute REML after chopping off N data points from the end of the data. Here is the 3dDeconvolve generated design matrix.

3dDeconvolve -input pb07.$subj.r*.scale+orig.HEAD                        \
    -censor censor_${subj}_combined_2.1D                                 \
    -ortvec mot_demean.r01.1D mot_demean_r01                             \
    -ortvec mot_demean.r02.1D mot_demean_r02                             \
    -ortvec mot_demean.r03.1D mot_demean_r03                             \
    -ortvec mot_demean.r04.1D mot_demean_r04                             \
    -ortvec mot_demean.r05.1D mot_demean_r05                             \
    -ortvec mot_deriv.r01.1D mot_deriv_r01                               \
    -ortvec mot_deriv.r02.1D mot_deriv_r02                               \
    -ortvec mot_deriv.r03.1D mot_deriv_r03                               \
    -ortvec mot_deriv.r04.1D mot_deriv_r04                               \
    -ortvec mot_deriv.r05.1D mot_deriv_r05                               \
    -polort 6                                                            \
    -num_stimts 3                                                        \
    -stim_times_AM1 1 stimuli/sub-xx_type1_dm.1D 'dmUBLOCK(-2)'    \
    -stim_label 1 type1                                                  \
    -stim_times_AM1 2 stimuli/sub-xx_type2_dm.1D 'dmUBLOCK(-2)'    \
    -stim_label 2 type2                                                  \
    -stim_times_AM1 3 stimuli/sub-xx_type3_dm.1D 'dmUBLOCK(-2)'    \
    -stim_label 3 type3                                                  \
    -jobs 16                                                             \
    -fout -tout -x1D X.xmat.1D -xjpeg X.jpg                              \
    -x1D_uncensored X.nocensor.xmat.1D                                   \
    -errts errts.${subj}                                                 \
    -x1D_stop                                                            \
    -bucket stats.$subj

Here is the 3dREMLfit code

3dREMLfit -matrix X.xmat.1D \
 -input "pb07.sub-xx.r01.scale+orig.HEAD pb07.sub-xx.r02.scale+orig.HEAD pb07.sub-xx.r03.scale+orig.HEAD pb07.sub-xx.r04.scale+orig.HEAD pb07.sub-xx.r05.scale+orig.HEAD" \
 -fout -tout -Rbuck stats.sub-xx_REML -Rvar stats.sub-xx_REMLvar \
 -Rerrts errts.sub-xx_REML -verb $*

To chop off some data, I could replace the above input by concatenating all runs and then taking whatever data point range.

3dTcat -prefix all_runs.$subj pb07.$subj.r*.scale+orig.HEAD

3dTcat -prefix all_runs.sub-epxxx.chunk1     \
    all_runs.sub-ep2169+orig.HEAD'[0..589]'

However, it is tricky to manipulate the design matrix 1D file, since 1) when the chopped length is longer than 1 run, there will be all 0 columns in the chopped design matrix 2) all 0 columns can be removed easily if all we need is the correct design matrix, but 3dREMLfit can not pass the 1D file header validation, such as ColumnLabels, GoodList, RunStart etc can not match with the modified design matrix.

Is there a convenient way to implement the above manipulation? Thanks a lot.

Best,
Zhengchen

2

Hi Zhengchen,

I would have to play around with that, too. But it seems much easier to rerun 3dDeconvolve -x1D_stop with truncated input, and have it generate the correct model. The data will need to be truncated anyway, to match the model.

Alternatively, you could cheat and modify just the -censor censor_${subj}_combined_2.1D option, zeroing out the last N rows. The betas should be the same, but I would have to ponder the statistics (or just ask @Gang! :). Of course, actually testing that would be an option, too, but then it might feel like actual work...

  • rick
1 Like
3

Hi Rick,

Thanks for your suggestion. I tried to hack the censor 1D file by setting the last N points to 0s. 3dDeconvolve indeed generates a different X.xmat.1D file, e.g., different ni_dimen, GoodList and the design matrix. However, I guess it is designed to remove just part of the data points, not a whole run. When N is larger than the length of one run, it still generates all zeros columns.

Probably, what I can do is when N < len(one run data point), use -censor; when N >= len(one run data point), first remove the run from 3dDeconvolve and then use -censor for the rest data points that require to be truncated. This can be done by a Python script to generate corresponding 3dDeconvolve and 3dREMLfit command.

Thanks,
Zhengchen

4

Actually, I was assuming it would be more complicated. Removing a final run should easy, just specify only the others with -input. Removing an earlier run is not so easy, since then your timing files would need changing.

The all zero columns are actually okay, though it means you would need an appropriate -GOFORIT option added.

  • rick
5

I only truncate from the end.

As far as I know, -GOFORIT just silences the warning message, but the actual design matrix will still contain all zero columns. I tried 3dDeconvolve -GOFORIT 99 and then 3dREMLfit, it failed with error ** ERROR: matrix column #28 is all zero!?

6

To be sure, what is wrong with having the all-zero columns? The purpose of that is just so the number of output volumes in the stats dataset is consistent, even with empty regressors. But the results from them should be all-zero, as well.

This seems okay.

  • rick
7

I guess when there is an all-zero column in the design matrix, GLM solver will face identifiability issues since the corresponding coefficient for the all-zero column can be any value (including 0). I don't know the code details, unless 3dREMLfit automatically ignores all-zero columns (solve GLM without all-zero columns and set 0 coefficients for the corresponding column name afterward), GLM with one or more empty columns would crush the model solving. Some algorithms may provide a solution until convergence, but it may not be stable. Anyway, I think 3dREMLfit indeed knows this issue, that's why it stopped with all-zero column errors.

8

It warns you of all-zero columns (under the expectation that it might be a mistake), bu if you run it (with -GOFORIT 99, say), the results for such volumes should be all zero as well, displaying that the program has no problem with it.

These programs should both support all-zero regressors without any trouble.

  • rick
9

Oh yes, I am sorry. I thought GOFORIT is only the option of 3dDeconvolve. I checked 3dREMLfit help, it also takes GOFORIT and explains how all-zero columns are handled. They are ignored and the coefficients are given as 0. I fit the model with GOFORIT, it works. Thanks a lot.

10

Great, thanks!

  • rick