Concerns about log-linear vs weighted log-linear vs direct monoexponential decay model fit

[dowdlelt]

Summary

From my understanding of the code, which may very well be wrong, we are performing a simple log-linear fit. While this is a very easy computation, I am concerned that it is somewhat incorrect, particularly in areas that drop out quickly. It may be that we are underestimating T2* and weighting the earlier echoes too heavily.

Additional Detail

For example, if you look at the thread here:
https://stackoverflow.com/questions/3433486/how-to-do-exponential-and-logarithmic-curve-fitting-in-python-i-found-only-poly
You note that fitting the exponential curve and using weights with log linear perform nearly identical, while log-linear fit with no weights comes up with a curve that drops very quickly.

I have done a bit of testing (have to get back to computer - still at scanner) that suggest this is true for our data, even in an eight echo test. While I'm working on getting those figures together, I wanted to start the discussion.

Note that this is seperate from the bi-exponential question in #212 , as the question in this issue is the accuracy of the (approximation) of the mono-exponential fit.

Next Steps

• Discuss - anyone else seen/thought about this?
• Add in a flag to perform slightly more computationally demanding curve/weighted fit

[dowdlelt]

It is also possible that it is estimating in the other direction, and overestimating T2*, leading to weighting late echoes too much. Need to rerun all the testing.

[rmarkello]

Hey @dowdlelt!

Sorry to just jump in out of the blue but I thought it might be relevant to link this back to up a previous post on the gitter by @handwerkerd. He references a similar concern about the T2* estimation and pointed out some work that @javiergcas had done in a separate repository implementing a curve fitting using scipy.optimize (see relevant codebits here). It definitely increases the computation time, so if using polyfit with a weighting parameter works then I think that might be a really good tradeoff!

[dowdlelt]

Thanks @rmarkello !

I knew this had been mentioned before! I searched through issues but forgot to go that far back in gitter. If I recall correctly, that was all related to attempting to use the goodness of fit throughout the subsequent steps.

In my brief testing, fitting a monoexponential does add time, and complexity but it wasn't a huge hit (for calculation on mean). The weighted fit is nearly indistinguishable from fitting scipy.optimize, while unweighted differs by up to 15ms on the T2* estimation.

Example image - this is the difference map between log-linear and weighted log-linear. This is weighted T2* - unweighted (Current method).

Range maxes out at 10 - so some areas in the orbitofrontal cortex are getting T2* estimates in the current method that are off by up to 15ms.

This is with an 8-echo data set. Performance is likely different (worse perhaps?) on a more typical range of echo times.

This is an example snippet of the code - its quite simple.

weights = np.sqrt(mean_array[:,vox_dim].flatten())
Mean_T2_w[vox_dim], Mean_S_weighted[vox_dim] = np.polyfit(TEs, log_mean_array[:,vox_dim].flatten(),1, w = weights)

where the weights are produced from the original array of the TEs, and polyfit is being fit to the log values, with order 1, and weights included.

Of course, the big question is if this actually is a meaningful difference.

[tsalo]

@dowdlelt Do you know if there's a basis for using sqrt(y) for the weights? Is it just what the original poster found to work well or is it based on the underlying math?

[dowdlelt]

@tsalo - I am not sure. I used it as it is what I found across a few places explaining weighted log-linear regression. On the actual data it also worked well and very closely matched the scipy.optimize results.

[jbteves]

Is mean_array squared earlier? If so maybe it's to guarantee all positive numbers.

[dowdlelt]

@jbteves negative - thats just the array size (n_echoes,n_voxels), that is just the mean of of each echoes timeseries.

[tsalo]

If it's something that came up in more than one place, then that's good enough for me. Do you have performance info? I.e., how much longer is the weighted version?

We currently have the argument --fitmode in the t2smap workflow, but not the tedana one, because we don't want to support volume-wise T2*/S0 estimation in the denoising workflow. We can add it into the tedana workflow with options of loglinear (default) and weighted-loglinear. We can also change the t2smap options from all and ts to loglinear (default), loglinear-volumewise (or something), and weighted-loglinear to make the two workflows consistent with one another.

[jbteves]

Alright, sorry I misunderstood. I'm pretty confused about what it could be then. Can you point me to one of the places you saw it explained?

[dowdlelt]

In my eight-echo dataset, performing the calculation on a mean volume, masked to just the brain -
log linear took 3.5 seconds
weighted took 9 seconds.
monoexponential optimize curve fit takes 38 seconds (log linear is also ran in this loop as well, to set initial guess for curve fit)

[tsalo]

You don't perform it on all of the data? We fit to all of the data in the log-linear fit currently implemented.

[dowdlelt]

@jbteves - I am only able to dig up the original post, though I thought I had found it in at least two different places - which gives me less confidence in it @tsalo .

@tsalo - At some point I fit it to all of the data, but I don't have estimates for how long that took because I was fitting all of these models, and a biexponential, and perhaps another log linear implementation. It took ages, as you might expect.

[dowdlelt]

@tsalo Found that code - running now, averaging about 7 to 8 timepoints per minute, consistent with estimated time above. So thats going to be much slower than the current tedana implementation. Not sure how this would change with just 3-echoes however.

It would be nice to run across multiple cores, but that's far over my head at the moment.

[jbteves]

Each timepoint is independent though, right? If so what you could do is multithread those jobs. The problem with the OpenMP stuff is that parallelizes with fairly reckless abandon and somehow worsens performance in our case. Users with a 4-8 core machine wouldn't see nearly as bad of a hit.

[dowdlelt]

oh yeah - its super easy in theory - its an independent job on each timepoint (or each voxel). However, I just managed to generate lots of errors on my first few tries and revealed I have no idea to multithread in python. (matlabs parfor was the full extent of my capabilities....)

[jbteves]

Multithreading is pretty tricky in most general-purpose languages. I fiddled with it in Python a while ago but it wasn't easy. Can you point me to a branch where you tried it or post a code snippet, and the errors? This issue is slowly centering around the idea that we should add that fit option, so I consider that useful info since somebody will probably open a PR eventually.

[dowdlelt]

I appreciate your high expectations - but I have no branches or code snippets. I hacked some things together in an ipython notebook, and didn't get anywhere with it and now it is gone.

@tsalo - total time for 220 timepoints, 8-echoes, weighted fit across 72013 voxels is 30 minutes. Thats a pretty substantial hit. In line with estimate, which suggest that my (poor) log linear implementation would take 10 minutes on this dataset (and scipy.optimze would be >2 hours).

[jbteves]

Yeah, but we can try to optimize around it. If the fit quality is better and we have concerns about it, this is something that only has to get done right once per dataset. It's not like ICA where there's a chance we fail and have to try again. Plus, if people can wait 6-8 hours for freesurfer I think waiting around an extra half hour for some good fits seems reasonable.
scipy.optimize does sound really, really slow though.

[tsalo]

I think I understand where we were getting our lines crossed. I don't mean fitting each volume independently (as in the FIT method), but rather fitting to all data, across time, simultaneously (as we do with numpy.linalg.lstsq in fit_decay).

To do that, we reshape the log data into (ExT, S) dimensions, and tile the (negative) echo times so that they're in (ExT,) shape. We then add a constant to the negative echo times so they're in (ExT, 2) shape. When we fit using least-squares, we get a constant and slope for each voxel, which we convert into T2* and S0.

My understanding is that polyfit will return the constant and slope by default (i.e., a and b from the SO example), so we shouldn't add a constant column to the negative echo times array. The same dimensions for the log data and negative echo times should be fine (except the negative echo times is 1D), and should give us the same results if we don't apply weighting as the standard OLS solution.

I also think it's weird to use the data for the weights, rather than the echo times. Is that really the way to go?

I'm thinking the approach should be to do this:

neg_tes = X[:, 1]  # ignore the old constants column
betas = np.polyfit(neg_tes, log_data, 1, w=(1 / np.sqrt(np.abs(neg_tes))))
t2s = 1. / betas[0, :].T # notice the switch in order
s0 = np.exp(betas[1, :]).T

We would put that here (or rather in the equivalent section of a new function):

tedana/tedana/decay.py

Lines 85 to 92 in 664b175

	log_data = np.log((np.abs(data[:, :echo + 1, :]) + 1).reshape(len(data), -1).T)
	# make IV matrix: intercept/TEs x (time series * echos)
	x = np.column_stack([np.ones(echo + 1), [-te for te in tes[:echo + 1]]])
	X = np.repeat(x, n_vols, axis=0)
	betas = np.linalg.lstsq(X, log_data, rcond=None)[0]
	t2s = 1. / betas[1, :].T
	s0 = np.exp(betas[0, :]).T

[dowdlelt]

I see, yes, I misunderstood. Removing the constant sounds right to me.

To me it makes sense to weight with signal intensities, since those change while echoes remain constant. Just a gut feeling though, and it matched monoexpo in testing

If I remember correctly the polyfit webpage has another example of weighting, which may be of interest.

[tsalo]

What shape are the weights in your implementation? (S,) or (S, E)?

The way I think about the problem is that we are trying to find T2* and S0 based on the data for a single voxel, and we have #TRs observations that we can use for that. So the data look like this figure:

I don't think it makes sense to weight observations based on the actual value, because that would mean downweighting certain time points. Instead, we want to downweight later echoes (or upweight earlier ones) when there is dropout and that slope flattens out, right? If that's the case, then we should probably use the echo time. Alternatively, we could use the value from the adaptive mask (number of good echoes) for this. We'd need to come up with a weighting scheme, but something like a 1 for good echoes and a 0.5 for bad ones could work.

[dowdlelt]

I'm not sure I follow - if we use the sqrt(signal_intensity) as the weight, this would give more weight to earlier echoes - as they have higher values. Wouldn't that satisfy your 'upweight earlier echoes' point?

Regarding your first question, weights are (8,1), as the sqrt of the mean signal at each echo time.

It may be that I'm viewing weighting as unrelated to the idea of precision, but rather related to creating a better 'monoexponential-like' fit. In this way, weighting based on values means we get better fits across the T2* across different voxels (which is very high) - whereas weighting based on precision means we may get better fits across the variability in T2*within a voxel over time, which is ~low (if it wasn't, BOLD would be easier!).

I don't feel that I am doing a very good job explaining this - as I just fiddled with my own implementation and you are working from the more reasonable point of the current code base.

[tsalo]

So you aren't weighting each voxel separately- I was definitely misunderstanding that, but it makes the code make sense again. We could hypothetically weight each echo in the array by the mean value of the data from that echo. It should definitely upweight earlier echoes.

But then, the weights would be for the whole brain, right? You end up with eight weights that are applied to all voxels simultaneously, rather than voxel-specific averages over time, correct?

[dowdlelt]

I believe I am still failing to explain myself well - perhaps I can clean up my notebook and share that. Is there a preferred way to share notebooks? It may take me a moment, but is perhaps better than propagating confusion.

[tsalo]

Sorry that I never responded to this. I think uploading your notebook to a fork is a good way to do it.

[dowdlelt]

I'm still getting settled in and figuring out things - hopefully I'll be able to get to this soon...

[dowdlelt]

I've implemented a direct curve fit (branch: curve_fit), in which t2smap (but not tedana) can take an additional flag:

                          dest='fittype',
                          action='store',
                          choices=['loglin', 'curvefit'],
                          help='Desired Fitting Method'
                                '"loglin" means that a linear model is fit'
                                ' to the log of the data, default'
                                '"curvefit" means that a more computationally'
                                'demanding monoexponential model is fit'
                                'to the raw data',
                          default='loglin')

Using curve fit means that the full fit for S0 (S0vG.nii) and T2star (t2svG) will reflect the monoexponential fit, rather than the log linear. The log linear fit is used to initialize the curvefit routine, and if fitting fails, the log lin value is used.

For an eight echo dataset this leads to changes that are sensible. Here is the S0 fit for the data with the log linear fit:

vs the S0 fit using curve fit:

To me, what stands out (beyond skull and such...) is the smoothness of intensity values in the brain. If S0 is to reflect the spin density of the voxel, then grey matter throughout the brain should be somewhat uniform. With the log linear fit, the S0 parameter is overestimated, as seen by the much higher values in low T2star regions such as the orbitofrontal cortex. Those regions don't have higher spin density, so they should match other regions. One bit of confusion here - this data likely had prescan norm on, as that is my default, so brightness due to coil sensitivity is not seen.

This of course has an effect on the T2star fits, left is the curvefir, right is log lin

Here the log lin as underestmated the T2 star value, and thus weighting is 'inaccurate'. Now, in practice this probably doesn't have a big effect on the timeseries (but I haven't tested that) but it could matter in regions prone to dropout.

Current implementation is very ugly, and doesn't respect the idea of fitting to differing numbers of echoes. But I offer it as a proof of principle. Because I never got around to uploading that jupyter notebook

[dowdlelt]

I will also take this opportunity to return to the idea of a weighted log linear fit - my previous method (which I encountered via google) gives values close to the curve fit method, but there is precedence for using weights calculated in a different manner, that I don't quite follow. Given here:

from here https://www.ncbi.nlm.nih.gov/pubmed/17374912

Per this link:
https://www.stat.cmu.edu/~cshalizi/350/lectures/18/lecture-18.pdf
It seems that weighting based on the inverse of the variance (per voxel) may be an approach.

The direct monoexponential curve fit isn't that much more computationally demanding unless a direct estimation of the S0 and T2star time series is desired, so may be preferable as a default

[tsalo]

@dowdlelt Do you have the weighted log-linear implemented too?

[dowdlelt]

sadly, no. I want to implement it as another fit option, as a sort of computational compromise - but I'm not sure that I can do it as the code does things now, I think it would instead be fit to a mean, rather than all data at once.

I'll try and put in in that branch in some manner, and post here when done.

[tsalo]

The curve-fitting method looks really cool. Do we have any way to validate it beyond interpreting the results qualitatively? Perhaps we could compare ME-EPI-based T2* estimates against a quantitative T2* map, if there are any datasets where both are available.

[dowdlelt]

I think there is a reason to have some trust in it - it is fitting the 'true' monoexponential to the raw data, rather than a log linear approximation.

Thats said - this is an excellent idea - I find the S0 result compelling, in the sense that there is no reason for the orbitofrontal cortex to have some excessively high spin density - but the effect on T2* is less - and it would be nice to have some ground truth.

I also need to do more testing with a fewer number of echoes - perhaps it is less exciting on more typical data....

[tsalo]

@dowdlelt, will you be able/willing to open a PR soon? I'd like to start trying this out on my own data. Plus, as @emdupre mentioned over in fMRIPrep, this might solve some of the issues with T2*-based coregistration.

[dowdlelt]

Yeah, I'll try and do that soon. I'll aim for today in fact, I think it is ready, in the sense that I'm not sure I can do much more with it without a lot of delay and work.

@tsalo - pull request started!