Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add AFNI 3dREMLfit for first-level estimation #171

Merged
merged 34 commits into from
Nov 19, 2020

Conversation

leej3
Copy link
Contributor

@leej3 leej3 commented Aug 30, 2019

Summary

Addresses #115

This is still a work in progress. I've created a FirstLevelModel nipype node that can be swapped in for the default fitlins one. To my knowledge, the estimator, which uses AFNI's 3dREMLfit program, can compute any F-tests and t-tests, with voxelwise REML estimation of the ARMA(1,1) temporal correlations (e.g., https://doi.org/10.1038/s41467-019-09230-w. However, there are some things that need to be considered for it to fully integrate with the fitlins workflow. Also, breaking it into separate nipype nodes and using the 3dREMLfit nipype interface rather than a subprocess call would probably be good.

For now. I have overwritten the fitlins.interfaces.nistats.FirstLevelModel _run_interface method and instead of using nistats.first_level_model.FirstLevelModel.fit I have written my own reml_fit as a method of fitlins.interfaces.afni.FirstLevelModel. These can be integrated as desired (I'm happy to help with some guidance) or persist as its own node/sub-workflow that can be swapped into the main fitlins workflow as requested by the user (likely by using the json rather than by the CLI change I’ve implemented here).

I overwrote fitlins.interfaces.utils.MergeAll because it checks for an equal number of contrast metadata elements and maps of each category. This allows me to make my way through a lot of the workflow but it just kicks the can down the road to the CollateWithMetadata class. The expectation that all elements have the same number is not met when one considers multirow F-tests (by multirow I am referring to F-tests specified with multiple condition weight vectors as part of the BIDS Stats Models specification json). I'm not sure what you want to do with this...

Other notes on the implementation

  • The main component of translation implemented here is to rewrite the design matrix and metadata from fitlins into a form that is consumable by 3dREMLfit. This form (described here, along with other useful features of the program) embeds the information in the NIML (XML-like) header of the AFNI design matrix file to guide the model fitting performed. It drops the need to specify many of the options of the 3dREMLfit CLI. I have hacked this together using string interpolation: it works but representing the metadata as XML and writing it to the file may be more elegant. I will look into that.

  • I'm not sure of the desired scope of this estimator. Would you prefer all results computed or a minimal output of the regression fit? As implemented, all statistical maps are computed using AFNI tools (3dREMLfit and then 3dPval). An alternative would be to tweak the data fed into 3dREMLfit (generate the 1d 'Y' object computed by nilearn) so that the output can be unmasked and the results maps can instead be computed using nistats.contrasts.compute_contrast. Regarding AFNI's computation of the maps...

    • I couldn't spot whether t-tests should be interpreted as 1- or 2-sided in the stat models specification (in the p-value computation). The AFNI team feels strongly about the fact that it should be 2-sided unless well-justified (and by far most FMRI-related cases are really 2-sided analyses)! I have used 3dPval to compute the p-value statistical maps, with the default being 2-sided. (For a 1-sided interpretation, these can be divided by 2)

    • z-score maps (the z statistic) are computed using a recent option added to 3dPval (afni/afni@7b41df7).

    • From the nistats source I think Bonferroni correction (for number of contrasts) is being applied during result map computation. This has not been added but can be within the reml_fit method in python if desired. Though perhaps this should not be a default and should instead be explicit in the stats model specification?

Possible additional features

  • adding multi-column stimuli (multiple basis functions). This requires only a small change in the code to implement. I will do that.
  • The minimize memory functionality within _run_interface is not implemented. I’m not sure what’s happening with that so would need help.
  • Similarly, caching of the fit function is not implemented.
  • Voxelwise covariates: at the run level, I believe AFNI would use the term voxelwise regressors. This 4D data orthogonalization could be implemented fairly easily with the 3dREMLfit “-dsort” option if required.
  • An Ordinary Least Squares could be computed in the same program run to compare the impact of the REML-estimated serial correlation pre-whitening. Or indeed the model could be fit quickly for a more basic analysis without ARMA, which is time-consuming (though better as documented in the help). Would this be a desired addition? It would likely provide useful comparisons with other non-REML estimators.
  • Many other tweaks to the computation that might be desirable. Perhaps these could be passed in as kwargs to reml_fit and injected into the command string/passed to the nipype 3d_REMLfit interface. Is this desirable or is the current constrained implementation and choices, or a variant thereof, acceptable?
  • For now, I’ve set the program to use a single thread to allow smooth operation of embarrassingly parallel computation of multiple BIDS runs on the same host. If this is used in a cluster setting that should probably be configurable to match resource specifications.
  • Are there other afni tools that would be especially useful for computing BIDS stats models? I'm happy to help out. Though perhaps waiting for the dust to settle after this addition would be good...
    In the stats models spec there is a description of how it attempts to expand model computation beyond the typical summary statistics paradigm. Are there currently implementations of this in fitlins or elsewhere? We (@afni-gangc) know of an approach that incorporates both effect estimate and the associated standard error as input for group analysis (3dMEMA in AFNI, FLAME is FSL, and also something in SPM along the same line).

Notes on neuroimaging processing

  • fitlins masks the data before fitting the regression. Typical AFNI processing does not do this. A basic extent mask is computed for the time series (a subset of all EPI voxels that have sufficient coverage over time to be able to meaningfully compute the regression and prevent issues during spatial blurring) but data outside the brain is not masked. As such, artifacts that are present and might be a cause for concern can be observed at the single level analysis output. Masking subsequently is applied before group level analyses/multiple comparisons corrections etc. With this said, if a mask file is supplied to the estimator, the 3dREMLfit calculations are done only on the voxels marked as being in the mask.
  • Also see the above note on t-test sidedness.
  • Censoring is an essential part of this sort of analysis. I had a brief chat with @effigies regarding its implementation in fitlins. Overall it seems like fitlins will be going with additional regressor columns in the design matrix to censor specific timepoints. Typically AFNI just removes the appropriate rows, but 3dREMLfit can process the time series dataset/design matrix either way. @afni-rwcox tested the results of both methods (column augmentation and row removal) and other than a few anomalous voxels (voxels identically zero except at the censored time points), the results were identical when fitting an ARMA model. One consideration for the choice would be that less RAM and CPU are required for the row removal method, because some volumes are excluded from the computation - this becomes significant when the fraction of censored time points goes over (say) 20% of the total.
  • I talked with @effigies about examples of reasonable applications of a multirow F-test. I got some useful ones:
    • modelling the hrf with multiple basis functions (TENT/CSPLIN functions in AFNI, FIR in SPM/FSL).
    • Main effects and interactions in an ANOVA design if the user wants to specify this
    • a quick and dirty goodness-of-fit test for the full model specification (e.g., do any of the specified conditions show a meaningful effect, worrying if the answer is "no" in the auditory cortex for an auditory task)

Queries about where the appropriate place for certain code-chunks.

  • Since AFNI volumewise statistical codes map onto nifti intent codes precisely, a recoder could be implemented in nibabel for this (as I have done here). Would this be a welcome pull request to nibabel?
  • On a more general note, would you accept pull requests for parsing other volumewise attributes to improve the API to BRIKs in nibabel? One bit of functionality that might be nice is for working with AFNI's "buckets". These are datasets that contain heterogeneous data volumes (e.g., beta coefficients and statistics for components in models). Splitting these into a list of homogenous volumes might be a useful method of the AFNIImage class. And implementing __iter__ for iterating through each element accessible by slicer[...,idx] might be nice...

@tyarkoni
Copy link
Collaborator

Holy shit, this fantastic, John! Had no idea you were working on this. Will read through the comments and try to answer any questions I can. Excited to use this! (@adelavega, feel free to experiment—it would be awesome to be able to give NeuroScout users the option to use AFNI for estimation).

@leej3
Copy link
Contributor Author

leej3 commented Aug 30, 2019 via email

@adelavega
Copy link
Collaborator

adelavega commented Aug 30, 2019

Very awesome. I guess in theory as far as NeuroScout is concerned the API should be plug and play interchangeable, correct? I'll take a closer look at the code next week!

Copy link
Collaborator

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for your patience. I've had a general look through, and made hopefully high level comments and questions. There are a lot of details to get to still, so I'll make more time on Monday.

fitlins/cli/run.py Outdated Show resolved Hide resolved
(24, "log10 p value", "Log10Pval"),
),
fields=("code", "label", "stat_code"),
)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pretty sure this is just nibabel.nifti1.intent_codes, restricted to the first 24?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The intent codes in nifti mirror the afni brik stat codes (they were copied from them). I don't know what's more intuitive to describe it as a stat code mapping or to add the brik labels as aliases to the niftit intent code mapping. I think moving this to nibabel would be good though.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since they're in NIFTI, I would probably just stick with that, as it's all official and whatnot. But if you wanted to add them to nibabel.brikhead, you're welcome to.

fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
nistats_flm.labels_, nistats_flm.results_, self.design_matrices_ = ([], [], [])
n_runs = len(run_imgs)
t0 = time.time()
for run_idx, run_img in enumerate(run_imgs):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it important to do these in series, or are they independent? If they're independent, we may be able to rewrite this to run just one, and use a MapNode to iterate.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

map node sounds good

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The workflow already uses boldfile as an iterfield in a map node. So the iteration is already done intelligently. the remlfit function (just like the nistats.fit function) has the ability to take a list of images though. I can remove that if you wish in the interest of clarity and brevity.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay. Let me have another read through before we refactor.

else:
params = [x for x in val[1].split(",")]
intent_info.append(
(STAT_CODES.label[val[0]], tuple([int(x) for x in params if x]))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
(STAT_CODES.label[val[0]], tuple([int(x) for x in params if x]))
(nb.nifti1.intent_codes.label[val[0]], tuple(int(x) for x in params if x))

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this will work. The AFNI stats code labels are very slightly different from the nifti intent codes.

fitlins/interfaces/afni.py Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
pval_bucket, pval_labels = load_bucket_by_prefix(tmpdir, pvals)
zscore_bucket, z_labels = load_bucket_by_prefix(tmpdir, zscores)
extra_bucket, extra_labels = load_bucket_by_prefix(tmpdir, extra_vars)
stat_bucket, stat_labels = load_bucket_by_prefix(tmpdir, glt)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

L332-L361: Is this stuff not already done in Remlfit? I assumed all of the parsing work should already be done. If it's not, it would make more sense to fix up that interface than to have a separate implementation here.

I guess we need to add a 3dPval interface. We can build that here or in nipype. (It'll get upstreamed eventually.)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good. That does seem like 30 lines of superfluous code. I didn't want to complicate things before a refactor. Would you have just defined and executed the single interface node within the class method?

I can take a hack at writing the 3dPval interface

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My thought here was to just use the interface, not wrapped in a Node(), which means it would just use the working directory of the Node that the whole interface is being run in.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great, will do.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have implemented this now. Let me know what you think

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps some of this code could be moved into the remlfit interface. The main problem is that remlfit is writing out images with different intent codes to the same file. These files need to be parsed and the results rewritten according to the categories required by FitLins.

And at some point the 3dPval interface can be moved to Nipype.

@effigies
Copy link
Collaborator

So I began to make code suggestions, but I think it makes more sense to step back and sketch out an overall API.

Your fitlins.interfaces.afni.FirstLevelModel is very nearly identical to fitlins.interfaces.nistats.FirstLevelModel up to and including the definition of a nistats.first_level_model.FirstLevelModel object. That says to me that what we really want there is to split that stuff out into something like a fitlins.interfaces.nistats.DesignMatrix interface that can be used directly.

Let's consider some set of abstract interfaces:

class DesignMatrixInterface(BaseInterface):
    class input_spec(BaseInterfaceInputSpec):
        session_info = traits.Dict(mandatory=True)
        bold_file = File(exists=True, mandatory=True)
    class output_spec(TraitedSpec):
        design_matrix = File()

class FirstLevelEstimatorInterface(BaseInterface):
    class input_spec(BaseInterfaceInputSpec):
        bold_file = File(exists=True, mandatory=True)
        mask_file = File(exists=True)
        smoothing_fwhm = traits.Float(desc='Full-width half max (FWHM) in mm for smoothing in mask')
        design_matrix = File(exists=True, mandatory=True)
        contrast_info = traits.List(traits.Dict)
    class output_spec(TraitedSpec):
        effect_maps = traits.List(File)
        variance_maps = traits.List(File)
        stat_maps = traits.List(File)
        zscore_maps = traits.List(File)
        pvalue_maps = traits.List(File)
        contrast_metadata = traits.List(traits.Dict)

Then we could simply say:

class nistats.DesignMatrix(DesignMatrixInterface):
    def _run_interface(self, runtime):
        import nibabel as nb
        from nistats import design_matrix as dm
        info = self.inputs.session_info
        img = nb.load(self.inputs.bold_file)
        vols = img.shape[3]

        if info['sparse'] not in (None, 'None'):
            sparse = pd.read_hdf(info['sparse'], key='sparse').rename(
                columns={'condition': 'trial_type',
                         'amplitude': 'modulation'})
            sparse = sparse.dropna(subset=['modulation'])  # Drop NAs
        else:
            sparse = None

        if info['dense'] not in (None, 'None'):
            dense = pd.read_hdf(info['dense'], key='dense')
            column_names = dense.columns.tolist()
            drift_model = None if 'cosine_00' in column_names else 'cosine'
        else:
            dense = None
            column_names = None
            drift_model = 'cosine'

        mat = dm.make_first_level_design_matrix(
            frame_times=np.arange(vols) * info['repetition_time'],
            events=sparse,
            add_regs=dense,
            add_reg_names=column_names,
            drift_model=drift_model,
        )

        mat.to_csv('design.tsv', sep='\t')
        self._results['design_matrix'] = os.path.join(runtime.cwd,
                                                      'design.tsv')
        return runtime


class nistats.FirstLevelEstimator(FirstLevelEstimatorInterface):
    def _run_interface(self, runtime):
        import nibabel as nb
        mat = pd.read_csv(self.inputs.design_matrix, sep='\t')
        img = nb.load(self.inputs.bold_file)
        mask_file = self.inputs.mask_file
        if not isdefined(mask_file):
            mask_file = None
        smoothing_fwhm = self.inputs.smoothing_fwhm
        if not isdefined(smoothing_fwhm):
            smoothing_fwhm = None
        flm = level1.FirstLevelModel(
            mask_img=mask_file, smoothing_fwhm=smoothing_fwhm)
        flm.fit(img, design_matrices=mat)

        effect_maps = []
        variance_maps = []
        stat_maps = []
        zscore_maps = []
        pvalue_maps = []
        contrast_metadata = []
        out_ents = self.inputs.contrast_info[0]['entities']
        fname_fmt = os.path.join(runtime.cwd, '{}_{}.nii.gz').format
        for name, weights, contrast_type in prepare_contrasts(
                self.inputs.contrast_info, mat.columns.tolist()):
            maps = flm.compute_contrast(weights, contrast_type, output_type='all')
            contrast_metadata.append(
                {'contrast': name,
                 'stat': contrast_type,
                 **out_ents}
                )

            for map_type, map_list in (('effect_size', effect_maps),
                                       ('effect_variance', variance_maps),
                                       ('z_score', zscore_maps),
                                       ('p_value', pvalue_maps),
                                       ('stat', stat_maps)):
                fname = fname_fmt(name, map_type)
                maps[map_type].to_filename(fname)
                map_list.append(fname)

        self._results['effect_maps'] = effect_maps
        self._results['variance_maps'] = variance_maps
        self._results['stat_maps'] = stat_maps
        self._results['zscore_maps'] = zscore_maps
        self._results['pvalue_maps'] = pvalue_maps
        self._results['contrast_metadata'] = contrast_metadata

        return runtime

And then your changes boil down roughly to your _reml_fit() function, which could now be afni.FirstLevelEstimator._run_interface.

I would recommend against using nistats.first_level_model.FirstLevelModel objects to store all of the design parameters. If those are things we need, I would rather make them an explicit part of the input specification, and pass them to FirstLevelEstimatorInterfaces.

Does that seem like a reasonable approach to you?

@leej3
Copy link
Contributor Author

leej3 commented Sep 27, 2019

Does that seem like a reasonable approach to you?

Yes, it looks good. Shall I make a first pass at it or do you want to do it and I can rebase and fix?

A problem I have noticed is that a variable number of output files are required. For example the number of stats_maps, effect_maps etc depend on the number of contrasts specified. This will get slightly more confusing when multi-row F tests (as defined above) are used because the count for each type of map written to disk will also differ.

I don't know whether the length of traited lists can be defined dynamically (based on parsing the model json). Or alternatively contrast could be an iterfield. Making it an iterfield would get the expected map count right (except for the case of multi-row F tests) but doesn't work with how the nistats prepare_contrasts function works. I'm not under the impression this would be better but thought I should flag it at this point rather than later...

@effigies
Copy link
Collaborator

I can do a first pass on Monday, but if you're planning on working the weekend, you're welcome to give it a shot.

@pep8speaks
Copy link

pep8speaks commented Oct 2, 2019

Hello @leej3, Thank you for updating!

Line 6:1: F401 'pandas as pd' imported but unused
Line 22:1: F401 '.utils.MergeAll' imported but unused
Line 68:9: F811 redefinition of unused 'pd' from line 6
Line 152:58: E241 multiple spaces after ','
Line 161:1: W293 blank line contains whitespace
Line 165:1: W293 blank line contains whitespace
Line 192:1: W293 blank line contains whitespace
Line 275:72: E202 whitespace before ')'
Line 275:100: E501 line too long (106 > 99 characters)
Line 302:100: E501 line too long (106 > 99 characters)
Line 302:107: W291 trailing whitespace
Line 305:36: W291 trailing whitespace
Line 306:100: E501 line too long (111 > 99 characters)
Line 312:32: E226 missing whitespace around arithmetic operator
Line 611:60: E231 missing whitespace after ','

Line 149:1: W293 blank line contains whitespace

Line 27:100: E501 line too long (147 > 99 characters)
Line 28:100: E501 line too long (143 > 99 characters)

Line 89:9: F811 redefinition of unused 'FirstLevelModel' from line 14
Line 128:100: E501 line too long (157 > 99 characters)
Line 234:17: E126 continuation line over-indented for hanging indent
Line 236:14: E131 continuation line unaligned for hanging indent
Line 237:14: E131 continuation line unaligned for hanging indent

To test for issues locally, pip install flake8 and then run flake8 fitlins.

Comment last updated at 2020-11-18 21:50:03 UTC

Copy link
Collaborator

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Had a couple minor wording/typo suggestions here, and I'm too lazy to cancel the review and make the main comment by itself.

fitlins/cli/run.py Outdated Show resolved Hide resolved
fitlins/cli/run.py Outdated Show resolved Hide resolved
from ..interfaces.visualizations import (
DesignPlot, DesignCorrelationPlot, ContrastMatrixPlot, GlassBrainPlot)
from ..interfaces.utils import MergeAll, CollateWithMetadata
from ..interfaces.afni import AFNIMergeAll
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are now being imported inside init_fitlins_wf. This works around some annoying new nipype behavior.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good to know

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was what was causing the packaging test failure.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Aha. That confused me.

@codecov-io
Copy link

codecov-io commented Oct 2, 2019

Codecov Report

Merging #171 (1b1fb35) into master (575bb58) will increase coverage by 2.73%.
The diff coverage is 86.24%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #171      +/-   ##
==========================================
+ Coverage   72.43%   75.16%   +2.73%     
==========================================
  Files          20       21       +1     
  Lines        1208     1506     +298     
  Branches      217      263      +46     
==========================================
+ Hits          875     1132     +257     
- Misses        234      261      +27     
- Partials       99      113      +14     
Flag Coverage Δ
afni 74.90% <84.89%> (?)
ds003 ?
nistats 58.49% <2.01%> (?)
pytest 12.61% <28.85%> (?)

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files Coverage Δ
fitlins/cli/run.py 82.83% <50.00%> (-1.02%) ⬇️
fitlins/interfaces/afni.py 86.45% <86.45%> (ø)
fitlins/interfaces/nistats.py 70.31% <100.00%> (+0.41%) ⬆️
fitlins/workflows/base.py 65.26% <100.00%> (+1.13%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 575bb58...1b1fb35. Read the comment docs.

@leej3
Copy link
Contributor Author

leej3 commented Oct 2, 2019

The nistats estimator automatically applies scaling. I have done the same albeit with a tweak to accommodate the 2d -> 4d switch in data shape. I don't know whether it would be preferable to specify this a transform in the model json instead.

I have assumed that the Tr (required for the afni design matrix) is in the design matrix index mat.index[1]. Is that a reasonable assumption?

@leej3
Copy link
Contributor Author

leej3 commented Oct 2, 2019

also I can't figure out why test_packaging fails

@effigies
Copy link
Collaborator

effigies commented Oct 2, 2019

I don't know

I have assumed that the Tr (required for the afni design matrix) is in the design matrix index mat.index[1]. Is that a reasonable assumption?

Assuming mat.index is in seconds not volumes, I'd think mat.index[1] - mat.index[0] would be more reliable, as there's always some kind of chance somebody's figured out how to start at t != 0. Do you know that it's in seconds? I haven't looked recently.

But if you have the nistats-built design-matrix, why do you need the TR again?

@leej3
Copy link
Contributor Author

leej3 commented Oct 2, 2019

Do you know that it's in seconds?

The tr in ds_00003 is 2.0, so I assume the unit is seconds.

But if you have the nistats-built design-matrix, why do you need the TR again?

I have to add it to the AFNI design matrix. session_info is an input to the design matrix node but not to the l1_model node. Rather than add it to the input spec I just wanted to extract it from that design matrix. I'll use the diff just in case.

so

mat = pd.read_csv(self.inputs.design_matrix, delimiter="\t", index_col=0)
t_r = mat.index[1] - mat.index[0]

@effigies
Copy link
Collaborator

effigies commented Oct 2, 2019

Ah, right, you need to translate to a new format. I think that's fine. The only potential issue is that there are experimental designs where there is no constant TR, e.g., clustered volume acquisition. We don't support this well yet, but it's a thing to keep in mind. Does AFNI handle situations like that?

@effigies effigies changed the title Add afni estimator Add AFNI 3dREMLfit for first-level estimation Nov 12, 2019
@effigies
Copy link
Collaborator

@leej3 Do we want to discuss this today? Sorry that it's taken so long to get around to it.

@leej3
Copy link
Contributor Author

leej3 commented Nov 14, 2019 via email

@effigies
Copy link
Collaborator

effigies commented Jan 6, 2020

Rerunning the failed test... Hopefully it was just a hiccup.

@effigies
Copy link
Collaborator

effigies commented Jan 8, 2020

Can you merge/rebase master?

@leej3
Copy link
Contributor Author

leej3 commented Jan 8, 2020 via email

@leej3
Copy link
Contributor Author

leej3 commented Jan 8, 2020

sorry I hadn't spotted your commit. You'll have to push that again.

@leej3
Copy link
Contributor Author

leej3 commented Nov 12, 2020

@effigies, are there any other changes you would like for this?

@effigies
Copy link
Collaborator

Thanks for this! I will look tomorrow.

@effigies
Copy link
Collaborator

I will look today.

Copy link
Collaborator

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Marking my progress. Meeting now, so will get back to it later. Feel free to respond now or wait for me to finish.

fitlins/cli/run.py Outdated Show resolved Hide resolved
fitlins/cli/run.py Outdated Show resolved Hide resolved
fitlins/cli/run.py Outdated Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
leej3 and others added 2 commits November 16, 2020 13:13
changes error ts arg/var name
explicitly uses runtime.cwd instead of tmpdir

Co-authored-by: Chris Markiewicz <effigies@gmail.com>
Copy link
Collaborator

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Small fixes here. I think this can go in and I'll make some issues for future improvements.

# modifying function as required and then append it to the
# appropriate output list type represented by map_list and write
# map_list to disk
for map_type, map_list, idx_list, mod_func in (
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mod_func doesn't appear to be used?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hmmm... that would be an error and mean that the effect variance is currently incorrectly the stdev. I'll check this.

fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
fitlins/interfaces/afni.py Outdated Show resolved Hide resolved
outmap.to_filename(fname)


class AFNIMergeAll(MergeAll):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A docstring explaining why this is happening would be helpful (I'm sure I knew, but now I don't), but this can be a new issue.

fitlins/interfaces/nistats.py Outdated Show resolved Hide resolved
fitlins/workflows/base.py Outdated Show resolved Hide resolved
fitlins/workflows/base.py Outdated Show resolved Hide resolved
adapt FirstLevelModel init for common api with AFNI estimator
tidy some code
improve nibabel usage

Co-authored-by: Chris Markiewicz <effigies@gmail.com>
@leej3
Copy link
Contributor Author

leej3 commented Nov 17, 2020

Ok thanks for taking a look. I'll take a pass through and push remaining suggestions:

  • Document MergeAll subclass
  • check that nibabel image api can be used (or switch to using dataobj)
  • propagate changes resulting from expectation that slicer will slice headers too.
  • fix or remove mod_func for results extraction (depending on whether np.square needs to be applied)

@leej3
Copy link
Contributor Author

leej3 commented Nov 17, 2020

@effigies

Document MergeAll subclass

So I admit the documentation is needed! I am pretty sure that my previous comment refers to this issue. In this subclass I overwrite the _list_outputs method to avoid the check of equal list lengths across fields. I don't have a sense of whether this violates the stats models IO spec or if there is an alternative solution that is more elegant. Do you have any suggestions?

check that nibabel image api can be used (or switch to using dataobj)

This looks fine. Nothing blew up.

propagate changes resulting from expectation that slicer will slice headers too.

This isn't required, the code handles either img or imgs fine. I'll leave get_afni_intent_info as is (and the full series is now passed to get_afni_intent_info_for_subvol as you suggested.

fix or remove mod_func for results extraction (depending on whether np.square needs to be applied)

@Shotgunosine and I spotted this a while back and then forgot to make the fix. He's working through that now.

@effigies
Copy link
Collaborator

effigies commented Nov 17, 2020

Document MergeAll subclass

So I admit the documentation is needed! I am pretty sure that my previous comment refers to this issue. In this subclass I overwrite the _list_outputs method to avoid the check of equal list lengths across fields. I don't have a sense of whether this violates the stats models IO spec or if there is an alternative solution that is more elegant. Do you have any suggestions?

I think let's just note in a new issue there's an inconsistent API. It may be that there's something in the nistats vs AFNI implementations, or it may be that the models you're using would require a similar change in nistats.

@leej3
Copy link
Contributor Author

leej3 commented Nov 17, 2020 via email

Copy link
Collaborator

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's do it. Thanks so much for pushing on this!

@effigies effigies merged commit 23c081a into poldracklab:master Nov 19, 2020
@tyarkoni
Copy link
Collaborator

Seconded—this is awesome!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

7 participants