BLR template configuration for a single-site model (Age, Age^2, Sex)

Hello,

I am currently setting up a voxel-wise normative model for diffusion MRI metrics across the adult lifespan. At the moment, my dataset comes from a single site (single batch).

My goal is to fit the following voxel-wise model: y ~β0 +β1*age +β2*age^2 +β3*sex

My input covariate matrix is structured as X = [age, sex] (with sex as a 0/1 dummy variable), and my batch effect array simply contains a single dummy batch for all subjects.

I would like to ask for a sanity check on my configuration to ensure it avoids redundant slope parameters for Age, while intentionally keeping the batch intercept available for future transfer learning.

Here is my current setup :

template_blr =  BLR(

    name = 'template_md_smooth', 

    basis_function_mean = PolynomialBasisFunction(basis_column = 0, degree  = 2) ,

    fixed_effect  = True,

    fixed_effect_slope  = False,

    fixed_effect_var_slope  = False,

    heteroskedastic =  True,

    warp_name   = None

)

 

# Normative Model

model_blr = NormativeModel(

template_regression_model = template_blr,

savemodel = True,

    evaluate_model = True,

    saveresults = True,

    saveplots = False,

    save_dir = output_blr,

    inscaler  = 'standardize',

    outscaler  = 'standardize',

    name = 'BLR_reference_model'

)

I set fixed_effect_slope = False to prevent the model from appending a duplicated ‘Age’ column in the design matrix, ensuring I get exactly one coefficient for the linear age effect, one for age^2, and one for sex, while I kept fixed_effect = True despite having a single site, to keep the batch intercept slot available for future transfer learning on ne datasets.

I set heteroskedastic = True to properly model the increasing variance expected in older individuals.

Any help or feedback would be greatly appreciated!

1 Like

Hi @elisabettaap and welcome to neurostars!

Have you considered using a generalized additive model (GAM) instead? It acts just like a GLM but allows for penalized splines that non-linearly model terms you choose (in your case, age). If longitudinal data, you can do a mixed-effect GAM (GAMM), or if cross-sectional, you may consider GAM for Location, Scale, and Shape (GAMLSS). GAMLSS tends to be the preferred method for making normative growth chart models and has been used in normative brain charting before: Brain charts for the human lifespan | Nature and White matter micro- and macrostructure brain charts for the human lifespan | Nature (second article is specific to dMRI).

For implementing these more complex statistical models efficiently in voxel-wise data, you can try the ModelArray/ModelArrayIO software (GitHub - PennLINC/ModelArrayIO: Companion converter software for ModelArray for converting data back and forth from the HDF5 file format. · GitHub and GitHub - PennLINC/ModelArray: ModelArray: an R package for statistical analysis of fixel-wise data and beyond · GitHub; and see the documentation linked in these pages).

Best,
Steven

Hi,

Thank you very much for your suggestion and references.

However, for this specific project, the pipeline is strictly built around the PCNtoolkit package and its Bayesian Linear Regression (BLR) implementation. Because of this, I really need to stick to the PCNtoolkit configuration.

If anyone familiar with this specific package could provide me feedback on what I mentioned in the original post, I would be truly grateful.

Best regards,
Elisabetta

Hi Elisabetta,

Nice to have you on Neurostars!

Your configuration looks pretty good to me. A few thoughts:

    • You might want to use a warped BLR (e.g. warp_name="WarpSinhArcsinh"), as dMRI metrics are usually skewed.

    • Keeping heteroskedastic=True for the reasons you mentioned is the right thing to do.

    • In pcntoolkit we usually choose categorical variables like sex as batch effects. You can then select fixed_effects=True and fixed_effect_var=True to model per sex mean + variance. With sex as batch effect PCNtoolkit can plot centiles of lets say total-gray-volume (y) to age (x) per sex and see how the y-x relationship changes for men vs women (a nice finding below is that men has a bit more gray volume than women):

      You can of course choose sex as a covariate. Then you get again per sex a mean and with heteroscedastic = True you model per sex a variance too.

    • fixed_effect=True: good idea to keep it on for future transfer and fixed_effect_slope is indeed redundant when you have a single dummy batch for all subjects.

    • Another nice idea: as dmri metrics usually are only positive, to avoid negative centiles going to values you can select y_transform="log"inside the NormativeModel.

Best,
Konstantinos