Move first-level stats up the chain with `block.get_contrasts()` about fitlins HOT 19 CLOSED

The CLI is usable. If you pip install -e ., you can run fitlins from the command line. Interactive usually goes something like:

$ ipython -i fitlins/base.py
In [1]: analysis = init('/data/bids/ds000114/model_noconfs.json', '/data/bids/ds000114', '/data/out/ds000114/derivatives/fmriprep')
In [2]: imgs = first_level(analysis, analysis.blocks[0], '/data/out/ds000114/derivatives/fitlins')
In [3]: mapping = second_level(analysis, analysis.blocks[1], '/data/out/ds000114/derivatives/fitlins')

from fitlins.

@tyarkoni If you've had any thoughts on what format/order you're expecting input files to be passed, that would be helpful. This matrix seems to be doing reordering, but if I wanted to use it to move from finger_vs_others to finger_vs_others.{test,retest}, I'd expect there to be some kind of sub-matrix that looks like:

                         finger_vs_others
finger_vs_others.test                   x
finger_vs_others.retest                 y

from fitlins.

I think your example has it the wrong way. Say you're at the subject level, and you're feeding in images from the run level. And say you have two runs, with all of the same images. Then the matrix returned by get_contrasts should look something like this:

condition               session	    finger_vs_others.test   finger_vs_others.retest
finger_vs_others        test        1.0                     0.0   
finger_vs_others        retest      0.0                     1.0   

The images coming in both have the same condition name, but the other entities (i.e., session) differ. It's the job of fitlins to make sure that the order of the images matches what you get from get_contrasts(). But once you've done that ordering, you should be able to just multiple the image data by the above matrix and get out images with the right names. I.e., something like:

dm = analysis.block['subject'].get_design_matrix().data
vox_data = ...  # Set up an n_voxels x n_images matrix
result = vox_data.dot(dm.values)  # this is now an n_voxels x n_contrasts matrix
for i, c in enumerate(dm.columns):
    img_data = result[:, i]
    # write out image to file, with condition name c
    ...

Am I missing a step? It seems to me that this should work all the way through, and I don't think you should need any custom code to handle the split transformation or any other transformation...

from fitlins.

I could use a get_contrasts matrix that looks like that. Right now it looks like what I posted above. I've resolved this for now by moving to a "mapping" dict that maps new column names onto old column names and relevant entities.

from fitlins.

Is your test workflow part of this code? If not, can you push a branch that has it? I just want to make sure I'm running the same example as you; will make it easier to debug.

from fitlins.

Actually, you can hold off. I figured out what I was missing. Should have a fix soon.

from fitlins.

Okay so I think the way this will work is that when you call get_contrasts(), you basically pass in (as the first argument) the results of a get_design_matrix() call--i.e., the full list of tuples. That will then set up the rows in the right order, and give you back a matrix in the format I specified. Does that work for you? Seems like the cleanest approach.

from fitlins.

I'm sorry, I missed the question at the end of that. I hope you haven't been waiting on me. That sounds like it should have enough information to make the transition between each level.

One question would be: does such a matrix lend itself to chaining in the sense that I concatenate (or similar) all of the contrast matrices in each level until I build one big second-level design matrix? If not, it's no big deal, but it would be nice to be able to do a two-shot model. I think that would resolve the NaN issue (nilearn/nistats#162) produced by saturated matrices at the intermediate level.

from fitlins.

No worries, I wasn't waiting on you--I have one other thing in the queue before getting to this, but it should happen tomorrow or Thursday.

Solving this problem in a sensible way required (or at least, induced) a fairly deep refactoring of the variables submodule. The upshot is that it will be easier to load variables for a single object (run, session, etc.), and yes, in practice, models will now be constructed/applied one at a time internally--so the way a higher-level model will be constructed is by concatenating lower-level ones. This sacrifices a little bit of efficiency (relative to the current approach of doing each level in one shot and then selecting out what one needs), but the gains in flexibility/sensibility are worth it, I think.

Will let you know once I push these changes.

from fitlins.

Just an update on my current thinking, in case you have some comments.

Regarding the contrast matrices, this is what I currently get:

In [1]: analysis = init('/data/bids/ds000114/model.json', '/data/bids/ds000114', '/data/out/ds000114/derivatives/fmriprep')
In [2]: block = analysis.blocks[1]
In [3]: cnames = [contrast['name'] for contrast in block.contrasts]
In [4]: block.get_contrasts(names=cnames)[0]
Out[4]: 
ContrastMatrixInfo(data=                         session_diff
finger_vs_others.test               1
finger_vs_others.retest            -1
finger_vs_others.test               1
finger_vs_others.retest            -1, index=  session subject            task
0  retest      01  fingerfootlips
1  retest      01  fingerfootlips
2    test      01  fingerfootlips
3    test      01  fingerfootlips, entities={'subject': '01', 'task': 'fingerfootlips'})

There are two issues here. The first is that there is a row that's indexed by retest that refers to test.finger_vs_others, and vice versa. The second is that there's still no clear relationship to the outputs of the previous stage (which were finger_vs_others). I need to know to parse the index, at which point we're back to the above issue where fitlins needs to be aware of all of the ways that transformed variable names may relate to their inputs. Instead, what would be ideal would be:

In [4]: block.get_contrasts(names=cnames)[0]
Out[4]: 
ContrastMatrixInfo(data=                  session_diff
finger_vs_others            -1
finger_vs_others             1, index=  session subject            task
0  retest      01  fingerfootlips
1    test      01  fingerfootlips, entities={'subject': '01', 'task': 'fingerfootlips'})

At which point I could do the following:

In [5]: data, idx, ents = block.get_contrasts(names=cnames)[0]
In [6]: for contrast_name, sub_ents in zip(data.index, idx.to_dict(orient='record')):
   ...:     layout.get(contrast=snake_to_camel(contrast_name), **sub_ents)

(This assumes that layout has been augmented to know of contrast entities.)

This seems to be the same issue from when we had a session-level block, where there was no way to identify finger_vs_others with test.finger_vs_others. I haven't digested the transforms code enough to have a feeling for how to get there, yet.

from fitlins.

I'm about to drive to Houston, so can't look into this at the moment. But does the approach I suggested on gutter for matching the inputs to the row names not work? I.e., look for a perfect match, and if it doesn't exist, split on '.' and match based on the entity index?

from fitlins.

The extra rows shouldn't be there though; maybe the transformations need to be moved level or something.

from fitlins.

Yeah, I can do that. Still feels kind of hacky, but okay. But it does give bad results with those extra rows.

Moving transformations back to their own "session" level again makes the "subject" level matrix look more sensible, but now the "session" level contrast matrix is huge and unwieldy, as there are now no actual contrasts to help constrain it to only variables of interest. For instance, every confound would get its own output.

from fitlins.

Hmm. Though maybe if we just treat the empty contrast list as a signal to do nothing, that might resolve the issue.

from fitlins.

Okay, so I was able to do this without maintaining my own mapping, using a layout of the fitlins derivatives and using the knowledge that the original contrast name must be one of a dot-separated series of names. It still feels hacky, but not less so than what I'd been doing.

I'm not using the "identity" contrasts at all. If you're not using them internally, or see some way they relate to eventually chaining, then I'm not sure they are used much.

I'm also not in any way using get_design_matrix() after the first level. Is there some way you anticipate that being used, or is that basically a first-level consideration?

Basically, I'm wondering if I'm failing to consider certain cases, and effectively coding too close to the prototype models.

from fitlins.

I think producing a contrast for every design matrix variable is really common; I almost always want to see those contrasts, even if they're not the focus of interest. So there should definitely be an option to add that automatically (and I think it will likely end up as part of the spec too). I think the more pertinent question is whether the identity_contrasts option should automatically include contrasts for all columns in the design matrix, or only for a subset (e.g., columns that are meant to be convolved with the HRF). I agree with you that confounds should almost certainly not be assigned to contrasts by default. What do you think about having the auto-contrast option only include variables named in hrf_variables in the JSON model spec? Or maybe it could be a string that takes values 'all' or 'hrf', with the latter being the default. The spec would have to be updated to match, so we would need to run this by other folks.

Re: post-first-level design matrices, that's a great question, and is something we've largely punted on so far when coming up with the draft BIDS-Model spec. There are definitely cases where people will need to do actual model fitting past the first level; for example, in the case of a whole-brain voxel-wise test for differences between two groups, but with continuous covariates added. That would be modeled at the dataset level, and the model would be something like y ~ group + 1|subject + covariate1 + covariate2 + ...

The thing that's kind of weird about this is that the current convention (which I implemented this way in the spec largely for consistency with other packages) is to basically embed the modeling into the contrast definition. I.e., instead of contrasts just being weighted sums of design matrix columns that create new columns (so, essentially, just another transformation), we use them to specify inferential tests. The more principled way to do it would be to separate the two. I.e., there would be nothing like a 'type' field in contrasts; instead, if you want to do a t-test on a particular condition, you'd have to specify that in the model section at the next level (and the actual specification would always be in terms of GLM/multiple regression; the t/F part would only come in when choosing what statistics to compute for the estimated model). So really, if we were being principled about it, I think a group-level test should look like this:

"models": [
  {
    "name": "group_diff.face",
    "variables": ["face", "age", "iq"],
    "statistics": ["t", "F", "R2"]
  },
  {
    "name": "group_diff.place",
    "variables": ["place", "age", "iq"],
    "statistics": ["t", "F", "R2"]
  },
  ...
]

The way we do it right now, we shove the stats stuff into contrasts, but that leaves us with no easy way to specify a model with subject-level covariates, as in the above case.

The upshot of all this is that maybe we need to rethink this part of the spec, and ultimately opt for the more verbose and redundant approach, even if it ends up being a bit of a pain in the ass for humans to read/write. Alternatively, we could expand the contrast definition to include a field like covariates, but I worry that that really starts to blur the line between what goes in model and what goes in contrasts.

Anyway, let me know what you think would make sense from a fitlins perspective.

from fitlins.

Oh, and to add to that: I still haven't added specification and handling of fixed and random effects to the spec, and that will definitely introduce cases where the design matrix will need to be manipulated post-run level. Ultimately, we will need to be able to get not one but two design matrices--one for the fixed effects and one for the random effects (where the latter is a binary matrix defining the grouping structure of observations). You definitely don't need to worry about this right now though, since there's no fMRI analysis software that can actually fit arbitrary mixed models properly. So it's a moot point implementation-wise, but I do think we need to ensure the API is compatible with planned future behavior, and get_design_matrix will eventually need to be called at other levels too.

from fitlins.

What do you think about having the auto-contrast option only include variables named in hrf_variables in the JSON model spec?

That seems reasonable. Unless we want to add "confounds" or "nuisance" or similar to the spec to indicate variables of non-interest. So HRF_variables are convolved, and set(variables) - set(nuisance) have contrasts generated. And the default contrast is the equivalent to FEAT's pe*.nii?

And pulling these individual variables up a level becomes a mean, e.g., across sessions?

from fitlins.

Anyway, let me know what you think would make sense from a fitlins perspective.

It looks like nistats is able to handle a design matrix in second-level models, but without an example model that requires using it, I don't know the bounds of its capabilities. So the fitlins perspective is somewhat fuzzy, on that matter.

It's not really clear to me how the statistic type makes sense as part of the model as opposed to the contrast. From a nistats perspective, this is a parameter for compute_contrast.

from fitlins.