enh - add Y transformation keyword argument in NormativeModel to enforce positivity

[contsili]

Closes #356

If Y<-1 then log(Y+1) cant be computed (ValueError) but I expect Y to always be >=-1 (since the functionality i added in this pr is only for data that are always positive). If we want to be conservative we can add though a check like:

# Validate that values are > -1 (required for log1p)
min_val = float(np.min(data[var].values))
if min_val <= -1.0:
    raise ValueError(
        f"y_transform='log1p' requires all {var} values > -1 but your {var} are not"
    )

We can still consider to add a test and/or an example.

[AuguB]

Will this solve the problem addressed in issue #365? As I see it, since the standardization happens after the log transform, this can and probably will still yield negative values, resulting in invalid input to WarpLog.

Additionally, why would we apply a log transform as a preparation for a log transform?

Some ideas:

• Explicitly disallowing standardize+WarpLog. A minmaxscaler can be used instead. Informative messages to the user can be helpful here.
• Adding a new scaler that only divides the data by its standard deviation, which amounts to a scaling around y=0; positive data stays positive. Such a scaler could be automatically applied when all of the following hold:

1. Y is strictly positive
2. BLR with Warplog as regressionmodel
3. The selected outscaler is standardize
   Then we can safely overwrite the outscaler with the new scaler (e.g. scale_only_standardize)

I do agree the option to add a log transform on Y as it is implemented now is a good idea, and it should be kept somehow. I just don't think it really solves the problem of using WarpLog + standardize. EDIT: I realize now that this PR is not intended to solve #365.

[contsili]

as you mentioned in the end of your previous comment, this pr is not meant to solve #365, but #356. #356 was proposed by @amarquand. It is meant to not allow for strictly positive response variables (like brain volumes or WM hypoinstensities) to have negative centiles (see for example, fig. 6(e) in [1]

For #365 for now we implemented something close to your first idea: see #366. However, your second idea is also a good one. We can make an issue to implement it.

[amarquand]

Indeed. We have a user that wants to apply the shash model to strictly positive data, and get the centiles back in the original space. But @AuguB is right that the interaction with the inscaler and outscaler needs to be carefully considered.

Has this been tested?

[contsili]

I wrote an example script which i have locally.

I used inscaler and outscaler = 'standardise', dataset fcon1000:

1. BLR(heteroskedastic=True)
Not enforced positivity:

Enforced positivity:

2. HBR() (default parameters)
Not enforced positivity

Enforced positivity:

@amarquand how do these results look?

[contsili]

I will add a test script for this feature and then we are ready to merge it

[contsili]

After discussion with @amarquand: I will implement two keyword arguments

1. y_transform = 'log'
2. y_transform = 'log1p'

In the case of 'log' we expect all centiles >=0 vs 'log1p' we expect centiles >=-1

[contsili]

from the tests the centile plots below are created (model parameters: BLR with linear basis, homoskedasticity, no warping) :
y-transform = None

y-transform = 'log'

y-transform = 'log1p'

We see that the fit is not good as we have non linear data with negative skewed and fit a very simple blr.

An interesting finding is that with the y-transform = 'log' all the centiles are positive but the 95th is quite far from the others. I belive that happens cause when the model’s predictions are made in log space and then convert them back to the original scale using exp(), the exp() "inflates" the differences between centiles.

The reason why the 95th centile is more inflated in the log case vs the log1p case is:
log1p(Y) maps this to [0, ...]
log(Y) maps this to [-inf, ...]

This wide range of log(Y) make the std of Y a lot higher than the log1p(Y), leading to an inflated centile = exp( Z * std + mean1)