HBR with SHASH likelihood

Welcome to this tutorial notebook that will go through the fitting and evaluation of Normative models with a Hierarchical Bayesian Regression (HBR) model using a SHASH likelihood.

Let’s jump right in.

Imports

import warnings
import logging


import pandas as pd
import matplotlib.pyplot as plt
from pcntoolkit import (
    HBR,
    BsplineBasisFunction,
    NormativeModel,
    NormData,
    load_fcon1000,
    SHASHbLikelihood,
    make_prior,
    plot_centiles_advanced,
    plot_qq,
    plot_ridge,
)

import numpy as np
import pcntoolkit.util.output
import seaborn as sns

sns.set_style("darkgrid")

# Suppress some annoying warnings and logs
pymc_logger = logging.getLogger("pymc")

pymc_logger.setLevel(logging.WARNING)
pymc_logger.propagate = False

warnings.simplefilter(action="ignore", category=FutureWarning)
pd.options.mode.chained_assignment = None  # default='warn'
pcntoolkit.util.output.Output.set_show_messages(False)

Load data

First we download a small example dataset from github.

# Download an example dataset
norm_data: NormData = load_fcon1000()

# Select only a few features
features_to_model = [
    "WM-hypointensities",
    "Right-Lateral-Ventricle",
    "Right-Amygdala",
    "CortexVol",
]
norm_data = norm_data.sel({"response_vars": features_to_model})

# Split into train and test sets
train, test = norm_data.train_test_split()

# Visualize the data
feature_to_plot = features_to_model[0]
df = train.to_dataframe()
fig, ax = plt.subplots(1, 2, figsize=(15, 5))

sns.countplot(data=df, y=("batch_effects", "site"), hue=("batch_effects", "sex"), ax=ax[0], orient="h")
ax[0].legend(title="Sex")
ax[0].set_title("Count of sites")
ax[0].set_xlabel("Site")
ax[0].set_ylabel("Count")


sns.scatterplot(
    data=df,
    x=("X", "age"),
    y=("Y", feature_to_plot),
    hue=("batch_effects", "site"),
    style=("batch_effects", "sex"),
    ax=ax[1],
)
ax[1].legend([], [])
ax[1].set_title(f"Scatter plot of age vs {feature_to_plot}")
ax[1].set_xlabel("Age")
ax[1].set_ylabel(feature_to_plot)

plt.show()

Creating a Normative model

A normative model has a regression model for each response variable. We provide a template regression model which is copied for each response variable.

A template regression model can be anything that extends the RegressionModel. We provide a number of built-in regression models, but you can also create your own.

Here we use the HBR class, which implements a Hierarchical Bayesian Regression model.

Likelihoods

HBR models are composed of a likelihood and a number of priors on the parameters of the likelihood. The PCNtoolkit offers a number of likelihood functions: 1. NormallLikelihood: Good for modeling data that is (approximately) normally distributed. 2. SHASHbLikelihood: Good for modeling data that is heavily skewed, or tailed. 3. BetaLikelihood: Good for modeling data that is bounded, e.g. between 0 and 1.

Likelihood parameters

Each of these likelihoods takes their own set of parameters, and for each, we have to set a prior: 1. NormalLikelihood: - mu: The mean of the normal distribution. - sigma: The standard deviation of the normal distribution. 2. SHASHbLikelihood: - mu: The mean of the skew-normal distribution. - sigma: The standard deviation of the skew-normal distribution. - epsilon: The skewness parameter of the skew-normal distribution. - delta: The tail thickness (or kurtosis) of the skew-normal distribution. 3. BetaLikelihood: - alpha: The shape parameter of the beta distribution. - beta: The scale parameter of the beta distribution.

Configuring likelihood parameters

Each likelihood parameter needs to be configured. The defaults should work reasonably well for most cases, at least where the data is standardized. Here’s a quick guide to configuring the likelihood parameters yourself, using the make_prior function.

1. Is your parameter a function of the covariates? If so, you set the linear parameter to True.

   1. If so, you can choose the basis expansion to use for the parameter: BSplineBasisFunction, LinearBasisFunction, or PolynomialBasisFunction.

   2. Also, determine whether the slope and intercept of the prior have a random effect or not. Here’s an example of a linear prior with a bspline basis expansion and a random effect in the intercept.

mu = make_prior('mu', linear=True, basis_function=BSplineBasisFunction(degree=3, nknots=5), intercept = make_prior('intercept_mu', random=True))

2. If your parameter is not a function of the covariates, you have to decide whether the parameter itself has a random effect or not. Here’s an example of a prior with a random effect.

epsilon = make_prior('epsilon', random=True)

3. Some parameters (such as sigma) need to be strictly positive, which we can enforce with a mapping. Here’s an example of a prior with a mapping to the positive real line.

# The mapping_params are (horizontal shift, scaling, vertical shift)
sigma = make_prior('sigma', mapping='softplus', mapping_params=(0, 5, 0))

# Mini demo of the mapping params
xsp = np.linspace(-7, 7, 100)
softplus = lambda x: np.log(1 + np.exp(x))
paramaterized_softplus = lambda x, a, b, c: softplus((x - a) / b) * b + c
plt.plot(xsp, paramaterized_softplus(xsp, 0, 1, 0), label="no mapping")
plt.plot(xsp, paramaterized_softplus(xsp, 1.5, 1, 0), label="horizontal shift of 1.5")
plt.plot(xsp, paramaterized_softplus(xsp, 0, 1, 1), label="vertical shift of 1")
plt.plot(xsp, paramaterized_softplus(xsp, 0, 2, 0), label="scale with a factor of 2")
plt.legend()
plt.show()

4. Any non-linear parameters can be further configured with dist_name and dist_params. Here’s an example of a prior with a gamma distribution.

alpha = make_prior('alpha', dist_name='gamma', dist_params=(1, 1))

We currently support the following distributions: - Normal - HalfNormal - LogNormal - Uniform - Gamma

The order of the parameters is important, and follows the order of the parameters in the corresponding distributions in PyMC.

Creating a HBR model

Here’s a thoroughly commented example of a HBR model with a SHASH Likelihood, which we will use to model our response variable.

# The SHASHb likelihood is a bit more flexible than the Normal likelihood, and takes four parameters, mu, sigma, epsilon, and delta.
# Mu and sigma fulfill the same role as in the Normal likelihood, namely the mean and standard deviation of the distribution.
# Epsilon and delta are parameters that control the skewness and kurtosis of the distribution.

# SHASHb model with fixed values for epsilon and delta

mu = make_prior(
    # Mu is linear because we want to allow the mean to vary as a function of the covariates.
    linear=True,
    # The slope coefficients are assumed to be normally distributed, with a mean of 0 and a standard deviation of 10.
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 10.0)),
    # The intercept is random, because we expect the intercept to vary between sites and sexes.
    intercept=make_prior(
        random=True,
        # Mu is the mean of the intercept, which is normally distributed with a mean of 0 and a standard deviation of 1.
        mu=make_prior(dist_name="Normal", dist_params=(0.0, 1.0)),
        # Sigma is the scale at which the intercepts vary. It is a positive parameter, so we have to map it to the positive domain.
        sigma=make_prior(dist_name="Normal", dist_params=(0.0, 1.0), mapping="softplus", mapping_params=(0.0, 3.0)),
    ),
    # We use a B-spline basis function to allow for non-linearity in the mean.
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
)
sigma = make_prior(
    # Sigma is also linear, because we want to allow the standard deviation to vary as a function of the covariates: heteroskedasticity.
    linear=True,
    # The slope coefficients are assumed to be normally distributed, with a mean of 0 and a standard deviation of 2.
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 2.0)),
    # The intercept is not random, because we assume the intercept of the variance to be the same for all sites and sexes.
    intercept=make_prior(dist_name="Normal", dist_params=(1.0, 1.0)),
    # We use a B-spline basis function to allow for non-linearity in the standard deviation.
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
    # We use a softplus mapping to ensure that sigma is strictly positive.
    mapping="softplus",
    # We scale the softplus mapping by a factor of 3, to avoid spikes in the resulting density.
    # The parameters (a, b, c) provided to a mapping f are used as: f_abc(x) = f((x - a) / b) * b + c
    # This basically provides an affine transformation of the softplus function.
    # a -> horizontal shift
    # b -> scaling
    # c -> vertical shift
    # You can leave c out, and it will default to 0.
    mapping_params=(0.0, 3.0),
)

epsilon = make_prior(
    # Epsilon is assumed to follow a normal distribution, with a mean of 0 and a standard deviation of 1.
    dist_name="Normal",
    dist_params=(0.0, 1.0),
)

delta = make_prior(
    # Delta is sampled from a normal distribution, with a mean of 1 and a standard deviation of 1, and then mapped to the positive real line using a softplus function.
    dist_name="Normal",
    dist_params=(1.0, 1.0),
    mapping="softplus",
    # We apply a softplus mapping to the delta parameter, to ensure that it is strictly positive.
    mapping_params=(
        0.0,  # Horizontal shift
        3.0,  # Scale for smoothness
        0.6,  # We need to provide a vertical shift as well, because the SHASH mapping goes a bit wild with low values for delta
    ),
)

shashb1_regression_model = HBR(
    name="template",
    cores=16,
    progressbar=True,
    draws=1500,
    tune=500,
    chains=4,
    nuts_sampler="nutpie",
    likelihood=SHASHbLikelihood(mu, sigma, epsilon, delta),
)

After specifying the regression model, we can configure a normative model.

A normative model has a number of configuration options: - savemodel: Whether to save the model after fitting. - evaluate_model: Whether to evaluate the model after fitting. - saveresults: Whether to save the results after evaluation. - saveplots: Whether to save the plots after fitting. - save_dir: The directory to save the model, results, and plots. - inscaler: The scaler to use for the input data. - outscaler: The scaler to use for the output data.

model = NormativeModel(
    # The regression model to use for the normative model.
    template_regression_model=shashb1_regression_model,
    # Whether to save the model after fitting.
    savemodel=True,
    # Whether to evaluate the model after fitting.
    evaluate_model=True,
    # Whether to save the results after evaluation.
    saveresults=True,
    # Whether to save the plots after fitting.
    saveplots=False,
    # The directory to save the model, results, and plots.
    save_dir="resources/hbr_SHASH/save_dir",
    # The scaler to use for the input data. Can be either one of "standardize", "minmax", "robminmax", "none"
    inscaler="standardize",
    # The scaler to use for the output data. Can be either one of "standardize", "minmax", "robminmax", "none"
    outscaler="standardize",
)

Fit the model

With all that configured, we can fit the model.

The fit_predict function will fit the model, evaluate it, and save the results and plots (if so configured).

After that, it will compute Z-scores and centiles for the test set.

All results can be found in the save directory.

model.fit_predict(train, test)

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 2 minutes

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.01	511
	2000	0	0.01	1023
	2000	0	0.02	1023
	2000	0	0.01	767

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 2 minutes

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.01	1023
	2000	0	0.01	767
	2000	0	0.01	511
	2000	0	0.01	1023

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 2 minutes

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.02	511
	2000	0	0.01	1023
	2000	0	0.01	255
	2000	0	0.01	1023

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 2 minutes

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.02	511
	2000	0	0.02	255
	2000	0	0.02	127
	2000	0	0.01	511

<xarray.NormData> Size: 98kB
Dimensions:            (observations: 216, response_vars: 4, covariates: 1,
                        batch_effect_dims: 2, centile: 5, statistic: 11)
Coordinates:
  * observations       (observations) int64 2kB 756 769 692 616 ... 751 470 1043
  * response_vars      (response_vars) <U23 368B 'WM-hypointensities' ... 'Co...
  * covariates         (covariates) <U3 12B 'age'
  * batch_effect_dims  (batch_effect_dims) <U4 32B 'sex' 'site'
  * centile            (centile) float64 40B 0.05 0.25 0.5 0.75 0.95
  * statistic          (statistic) <U8 352B 'EXPV' 'MACE' ... 'SMSE' 'ShapiroW'
Data variables:
    subject_ids        (observations) object 2kB 'Munchen_sub96752' ... 'Quee...
    Y                  (observations, response_vars) float64 7kB 2.721e+03 .....
    X                  (observations, covariates) float64 2kB 63.0 ... 23.0
    batch_effects      (observations, batch_effect_dims) <U17 29kB 'F' ... 'Q...
    Z                  (observations, response_vars) float64 7kB 0.5221 ... -...
    centiles           (centile, observations, response_vars) float64 35kB 1....
    logp               (observations, response_vars) float64 7kB -1.623 ... -...
    Yhat               (observations, response_vars) float64 7kB 2.362e+03 .....
    statistics         (response_vars, statistic) float64 352B 0.3927 ... 0.995
Attributes:
    real_ids:                       True
    is_scaled:                      False
    name:                           fcon1000_test
    unique_batch_effects:           {np.str_('sex'): [np.str_('F'), np.str_('...
    batch_effect_counts:            defaultdict(<function NormData.register_b...
    covariate_ranges:               {np.str_('age'): {'mean': np.float64(28.2...
    batch_effect_covariate_ranges:  {np.str_('sex'): {np.str_('F'): {np.str_(...

xarray.NormData

• Dimensions:
  • observations: 216
  • response_vars: 4
  • covariates: 1
  • batch_effect_dims: 2
  • centile: 5
  • statistic: 11
• Coordinates: (6)
  • observations
    (observations)
    int64
    756 769 692 616 ... 751 470 1043

    array([ 756,  769,  692, ...,  751,  470, 1043], shape=(216,))

  • response_vars
    (response_vars)
    <U23
    'WM-hypointensities' ... 'Cortex...

    array(['WM-hypointensities', 'Right-Lateral-Ventricle', 'Right-Amygdala',
           'CortexVol'], dtype='<U23')

  • covariates
    (covariates)
    <U3
    'age'

    array(['age'], dtype='<U3')

  • batch_effect_dims
    (batch_effect_dims)
    <U4
    'sex' 'site'

    array(['sex', 'site'], dtype='<U4')

  • centile
    (centile)
    float64
    0.05 0.25 0.5 0.75 0.95

    array([0.05, 0.25, 0.5 , 0.75, 0.95])

  • statistic
    (statistic)
    <U8
    'EXPV' 'MACE' ... 'SMSE' 'ShapiroW'

    array(['EXPV', 'MACE', 'MAPE', 'MSLL', 'NLL', 'R2', 'RMSE', 'Rho', 'Rho_p',
           'SMSE', 'ShapiroW'], dtype='<U8')

• Data variables: (9)
  • subject_ids
    (observations)
    object
    'Munchen_sub96752' ... 'Queensla...

    array(['Munchen_sub96752', 'NewYork_a_sub18638', 'Leiden_2200_sub87320',
           'ICBM_sub47658', 'AnnArbor_b_sub45569', 'Beijing_Zang_sub18960',
           'Leiden_2200_sub18456', 'Berlin_Margulies_sub27711',
           'Beijing_Zang_sub87776', 'Milwaukee_b_sub63196',
           'Beijing_Zang_sub07144', 'Atlanta_sub76280',
           'Beijing_Zang_sub40037', 'Cambridge_Buckner_sub17737',
           'ICBM_sub89049', 'ICBM_sub55656', 'Oulu_sub45566',
           'Beijing_Zang_sub89088', 'Atlanta_sub16563',
           'Cambridge_Buckner_sub51172', 'Oulu_sub98739',
           'Queensland_sub49845', 'Cambridge_Buckner_sub84256',
           'Cleveland_sub80263', 'ICBM_sub16607', 'Newark_sub46570',
           'NewYork_a_sub88286', 'Cambridge_Buckner_sub02591',
           'Oulu_sub66467', 'Beijing_Zang_sub74386', 'Newark_sub55760',
           'ICBM_sub30623', 'Oulu_sub68752', 'Leiden_2180_sub19281',
           'Beijing_Zang_sub50972', 'Beijing_Zang_sub85030',
           'Milwaukee_b_sub36386', 'Baltimore_sub31837', 'PaloAlto_sub84978',
           'Oulu_sub01077', 'NewYork_a_ADHD_sub54828', 'PaloAlto_sub96705',
           'Cambridge_Buckner_sub40635', 'ICBM_sub66794',
           'Beijing_Zang_sub46541', 'Beijing_Zang_sub87089',
           'Pittsburgh_sub97823', 'Beijing_Zang_sub98617', 'ICBM_sub92028',
    ...
           'Leiden_2200_sub04484', 'Beijing_Zang_sub80163', 'ICBM_sub02382',
           'Cambridge_Buckner_sub77435', 'NewYork_a_sub54887',
           'Oulu_sub85532', 'Baltimore_sub73823', 'Beijing_Zang_sub29590',
           'Oulu_sub99718', 'Beijing_Zang_sub08455', 'Beijing_Zang_sub85543',
           'Cambridge_Buckner_sub45354', 'Beijing_Zang_sub07717',
           'Baltimore_sub76160', 'Beijing_Zang_sub17093',
           'AnnArbor_b_sub90127', 'SaintLouis_sub73002',
           'Queensland_sub93238', 'Cleveland_sub34189',
           'Cambridge_Buckner_sub89107', 'Atlanta_sub75153',
           'NewYork_a_ADHD_sub73035', 'Cambridge_Buckner_sub59434',
           'Milwaukee_b_sub44912', 'Cleveland_sub46739', 'Oulu_sub20495',
           'SaintLouis_sub28304', 'Cambridge_Buckner_sub35430',
           'Oulu_sub86362', 'Newark_sub58526', 'Leiden_2180_sub12255',
           'ICBM_sub48210', 'Cambridge_Buckner_sub77989',
           'Berlin_Margulies_sub75506', 'NewYork_a_sub29216',
           'Beijing_Zang_sub05267', 'AnnArbor_b_sub18546', 'Oulu_sub75620',
           'AnnArbor_b_sub30250', 'Berlin_Margulies_sub86111',
           'Beijing_Zang_sub89592', 'Beijing_Zang_sub68012',
           'NewYork_a_sub50559', 'Munchen_sub66933',
           'Cambridge_Buckner_sub59729', 'Queensland_sub86245'], dtype=object)

  • Y
    (observations, response_vars)
    float64
    2.721e+03 1.289e+04 ... 5.035e+05

    array([[2.72140000e+03, 1.28916000e+04, 1.43940000e+03, 4.57858328e+05],
           [1.14310000e+03, 9.91910000e+03, 1.64970000e+03, 5.26780362e+05],
           [9.55800000e+02, 7.47730000e+03, 1.83850000e+03, 4.95744471e+05],
           [1.47390000e+03, 1.43021000e+04, 1.86770000e+03, 5.85303839e+05],
           [7.57800000e+02, 4.11930000e+03, 1.32500000e+03, 3.33111552e+05],
           [8.71100000e+02, 5.03090000e+03, 1.90730000e+03, 5.10794940e+05],
           [1.20730000e+03, 1.78664000e+04, 2.02220000e+03, 5.50533325e+05],
           [5.95000000e+02, 5.00790000e+03, 2.01070000e+03, 4.67673977e+05],
           [6.82400000e+02, 7.28660000e+03, 1.45630000e+03, 4.60129533e+05],
           [4.45100000e+02, 5.74290000e+03, 1.47450000e+03, 4.44494817e+05],
           [1.62000000e+03, 3.71370000e+03, 2.00110000e+03, 5.59424624e+05],
           [6.02800000e+02, 5.30120000e+03, 1.36100000e+03, 4.21551234e+05],
           [1.43250000e+03, 4.42970000e+03, 1.65080000e+03, 5.19842763e+05],
           [1.90820000e+03, 3.57810000e+03, 1.88370000e+03, 5.06679262e+05],
           [1.83400000e+03, 3.27190000e+03, 2.05120000e+03, 5.35569987e+05],
           [4.59600000e+02, 3.98580000e+03, 1.45470000e+03, 4.67607555e+05],
           [1.21000000e+03, 8.72130000e+03, 1.71430000e+03, 5.30904612e+05],
           [8.45900000e+02, 6.59310000e+03, 1.61900000e+03, 5.09371867e+05],
           [9.95200000e+02, 7.04020000e+03, 1.99490000e+03, 4.60068379e+05],
           [1.73470000e+03, 4.01480000e+03, 1.51620000e+03, 4.87269373e+05],
    ...
           [7.85800000e+02, 5.70900000e+03, 1.47480000e+03, 4.53982166e+05],
           [2.24010000e+03, 4.36660000e+03, 2.04210000e+03, 5.58453123e+05],
           [7.58100000e+02, 6.52980000e+03, 1.56730000e+03, 4.73575183e+05],
           [1.44050000e+03, 6.70530000e+03, 1.20540000e+03, 3.82788491e+05],
           [8.18600000e+02, 9.38330000e+03, 1.96740000e+03, 5.02713911e+05],
           [3.76990000e+03, 1.58644000e+04, 1.79170000e+03, 5.12490348e+05],
           [8.80200000e+02, 4.37020000e+03, 1.75520000e+03, 4.37300069e+05],
           [8.23900000e+02, 6.37900000e+03, 1.57650000e+03, 5.67331908e+05],
           [2.11390000e+03, 1.07225000e+04, 1.84380000e+03, 5.12273764e+05],
           [7.41900000e+02, 8.80170000e+03, 1.60640000e+03, 4.91973562e+05],
           [1.33390000e+03, 6.98000000e+03, 1.74850000e+03, 4.78907154e+05],
           [7.07300000e+02, 5.68070000e+03, 1.53450000e+03, 4.74077083e+05],
           [1.13410000e+03, 5.59220000e+03, 1.62620000e+03, 4.54163909e+05],
           [4.38600000e+02, 6.33000000e+03, 1.59670000e+03, 4.68067037e+05],
           [9.66300000e+02, 9.21550000e+03, 1.78250000e+03, 5.09199708e+05],
           [4.24300000e+02, 4.51110000e+03, 1.70200000e+03, 5.26635258e+05],
           [6.04700000e+02, 7.59080000e+03, 1.69930000e+03, 5.20499663e+05],
           [2.34320000e+03, 1.71923000e+04, 1.79380000e+03, 4.86680791e+05],
           [2.72170000e+03, 6.08600000e+03, 2.32470000e+03, 6.10402006e+05],
           [7.03500000e+02, 1.07003000e+04, 1.67620000e+03, 5.03535771e+05]])

  • X
    (observations, covariates)
    float64
    63.0 23.27 22.0 ... 72.0 23.0 23.0

    array([[63.  ],
           [23.27],
           [22.  ],
           [42.  ],
           [63.  ],
           [23.  ],
           [21.  ],
           [26.  ],
           [21.  ],
           [49.  ],
           [20.  ],
           [23.  ],
           [20.  ],
           [26.  ],
           [35.  ],
           [21.  ],
           [22.  ],
           [19.  ],
           [34.  ],
           [18.  ],
    ...
           [21.  ],
           [20.  ],
           [22.  ],
           [25.  ],
           [25.  ],
           [73.  ],
           [22.  ],
           [28.  ],
           [29.06],
           [19.  ],
           [20.  ],
           [22.  ],
           [19.  ],
           [24.  ],
           [21.  ],
           [24.  ],
           [22.79],
           [72.  ],
           [23.  ],
           [23.  ]])

  • batch_effects
    (observations, batch_effect_dims)
    <U17
    'F' 'Munchen' ... 'M' 'Queensland'

    array([['F', 'Munchen'],
           ['M', 'NewYork_a'],
           ['F', 'Leiden_2200'],
           ['M', 'ICBM'],
           ['F', 'AnnArbor_b'],
           ['M', 'Beijing_Zang'],
           ['M', 'Leiden_2200'],
           ['F', 'Berlin_Margulies'],
           ['F', 'Beijing_Zang'],
           ['F', 'Milwaukee_b'],
           ['M', 'Beijing_Zang'],
           ['F', 'Atlanta'],
           ['F', 'Beijing_Zang'],
           ['F', 'Cambridge_Buckner'],
           ['M', 'ICBM'],
           ['F', 'ICBM'],
           ['M', 'Oulu'],
           ['F', 'Beijing_Zang'],
           ['M', 'Atlanta'],
           ['F', 'Cambridge_Buckner'],
    ...
           ['F', 'SaintLouis'],
           ['M', 'Cambridge_Buckner'],
           ['F', 'Oulu'],
           ['F', 'Newark'],
           ['M', 'Leiden_2180'],
           ['M', 'ICBM'],
           ['F', 'Cambridge_Buckner'],
           ['M', 'Berlin_Margulies'],
           ['M', 'NewYork_a'],
           ['F', 'Beijing_Zang'],
           ['M', 'AnnArbor_b'],
           ['F', 'Oulu'],
           ['F', 'AnnArbor_b'],
           ['F', 'Berlin_Margulies'],
           ['M', 'Beijing_Zang'],
           ['F', 'Beijing_Zang'],
           ['M', 'NewYork_a'],
           ['M', 'Munchen'],
           ['M', 'Cambridge_Buckner'],
           ['M', 'Queensland']], dtype='<U17')

  • Z
    (observations, response_vars)
    float64
    0.5221 0.3941 ... -1.314 -1.141

    array([[ 5.22111034e-01,  3.94099738e-01, -2.66792814e-01,
             4.06460635e-01],
           [ 4.22466208e-01,  9.65448035e-01, -1.23115188e+00,
            -2.75083913e-01],
           [ 3.18291793e-01,  4.04732097e-01,  5.45620314e-01,
             2.88865481e-01],
           [ 5.69111007e-01,  1.36465454e+00, -2.58229393e-01,
             1.48401663e+00],
           [-2.02047143e+00, -1.70824599e+00, -2.79573577e-01,
            -1.90516036e+00],
           [-5.72758421e-01, -5.16773934e-01,  3.06331316e-01,
            -5.25596222e-01],
           [ 6.75159082e-01,  2.69725411e+00,  4.09408581e-01,
             3.40759379e-01],
           [-2.58377585e-02, -7.19266423e-01,  1.68705934e+00,
            -6.77610473e-01],
           [-1.09910335e+00,  6.32661043e-01, -9.32247384e-01,
            -6.61896768e-01],
           [-1.64782803e+00, -3.09579804e-01, -6.85539790e-01,
             2.02469096e-01],
    ...
           [ 3.23383225e-02,  4.56707221e-01,  3.94488157e-01,
            -4.75028151e-01],
           [-2.28010369e+00,  1.25793914e-01, -3.72871224e-01,
            -5.86864722e-01],
           [-3.68727350e-01,  2.77777856e-01,  8.28827400e-01,
             1.39496191e-01],
           [-5.46991665e-01, -3.98633468e-02, -1.89215774e-01,
            -7.79284242e-01],
           [-2.02910739e-01,  9.55953072e-01, -2.96613602e-01,
            -7.00202583e-01],
           [-2.75327977e+00, -4.66405172e-01,  3.47307354e-01,
             1.29650564e+00],
           [-1.66030307e+00,  3.97481554e-01, -9.24188289e-01,
            -4.76654211e-01],
           [-3.08981685e-01,  5.44325984e-01,  9.89606313e-01,
             2.91387578e-01],
           [ 2.46018649e+00,  2.29349516e-03,  1.98860470e+00,
             3.02005316e+00],
           [-9.65105272e-01,  1.16469286e+00, -1.31411301e+00,
            -1.14085260e+00]])

  • centiles
    (centile, observations, response_vars)
    float64
    1.004e+03 4.959e+03 ... 6.059e+05

    array([[[1.00398750e+03, 4.95918729e+03, 1.23858886e+03, 3.86736055e+05],
            [6.08602095e+02, 3.81136404e+03, 1.58699256e+03, 4.79606192e+05],
            [4.61188041e+02, 3.73382051e+03, 1.44083473e+03, 4.27724793e+05],
            ...,
            [1.09822390e+03, 5.73774119e+03, 1.30038803e+03, 4.16559291e+05],
            [9.34238404e+02, 3.51395053e+03, 1.60910669e+03, 4.41917504e+05],
            [5.79942119e+02, 3.76219376e+03, 1.62513992e+03, 4.86708761e+05]],
    
           [[1.52906096e+03, 7.38744766e+03, 1.37175546e+03, 4.18879836e+05],
            [7.86085048e+02, 4.98137572e+03, 1.74054148e+03, 5.12383721e+05],
            [6.40568118e+02, 4.86249913e+03, 1.59339913e+03, 4.60415329e+05],
            ...,
            [1.91753230e+03, 8.82741884e+03, 1.45358310e+03, 4.50275711e+05],
            [1.11202673e+03, 4.67512909e+03, 1.76241305e+03, 4.74666069e+05],
            [7.57730442e+02, 4.92337233e+03, 1.77844629e+03, 5.19457326e+05]],
    
           [[2.07580988e+03, 1.03198971e+04, 1.48322514e+03, 4.42912966e+05],
            [9.70933244e+02, 6.39440744e+03, 1.86911066e+03, 5.36897209e+05],
            [8.27397246e+02, 6.22553695e+03, 1.72114396e+03, 4.84863169e+05],
            ...,
            [2.77052962e+03, 1.25581568e+04, 1.58180962e+03, 4.75490524e+05],
            [1.29719402e+03, 6.07747726e+03, 1.89077914e+03, 4.99157770e+05],
            [9.42897736e+02, 6.32572049e+03, 1.90681237e+03, 5.43949027e+05]],
    
           [[2.94786482e+03, 1.51168517e+04, 1.61037303e+03, 4.67976306e+05],
            [1.26591979e+03, 8.70543430e+03, 2.01572517e+03, 5.62464780e+05],
            [1.12554392e+03, 8.45480271e+03, 1.86681796e+03, 5.10361722e+05],
            ...,
            [4.13092024e+03, 1.86606901e+04, 1.72803919e+03, 5.01790603e+05],
            [1.59268967e+03, 8.37103215e+03, 2.03716196e+03, 5.24702508e+05],
            [1.23839338e+03, 8.61927538e+03, 2.05319520e+03, 5.69493765e+05]],
    
           [[4.65019278e+03, 2.41226868e+04, 1.79662568e+03, 5.03652340e+05],
            [1.84205824e+03, 1.30430486e+04, 2.23036251e+03, 5.98856188e+05],
            [1.70784407e+03, 1.26390713e+04, 2.08007745e+03, 5.46654607e+05],
            ...,
            [6.78643678e+03, 3.01176185e+04, 1.94218858e+03, 5.39228041e+05],
            [2.16982050e+03, 1.26758769e+04, 2.25145994e+03, 5.61061376e+05],
            [1.81552422e+03, 1.29241201e+04, 2.26749317e+03, 6.05852634e+05]]],
          shape=(5, 216, 4))

  • logp
    (observations, response_vars)
    float64
    -1.623 -1.648 ... -1.373 -1.177

    array([[-1.62282548e+00, -1.64785626e+00, -6.29004113e-01,
            -7.04684293e-01],
           [-4.37893961e-01, -1.50704043e+00, -1.24324841e+00,
            -6.16205101e-01],
           [-3.53409858e-01, -8.94472096e-01, -9.82346850e-01,
            -6.67329860e-01],
           [-9.91335914e-01, -2.35266878e+00, -7.68952127e-01,
            -1.77734024e+00],
           [-2.46130899e+00, -1.80307654e+00, -6.15292810e-01,
            -2.29610528e+00],
           [ 1.69843795e-01, -2.76843190e-01, -8.49856792e-01,
            -6.82519637e-01],
           [-7.19451109e-01, -4.77541627e+00, -9.05046781e-01,
            -6.85689292e-01],
           [-6.04283284e-02, -3.85498427e-01, -2.26037994e+00,
            -8.00264151e-01],
           [-5.70768126e-02, -1.08106993e+00, -9.50155776e-01,
            -7.47096740e-01],
           [-1.34655206e+00, -7.91392769e-01, -8.05409464e-01,
            -6.93195709e-01],
    ...
           [-1.44179101e-01, -8.85692449e-01, -8.92070483e-01,
            -6.88624281e-01],
           [-1.90784482e+00, -6.47669808e-01, -7.26611956e-01,
            -7.10456058e-01],
           [ 4.77865405e-02, -6.97197063e-01, -1.18036135e+00,
            -6.40334051e-01],
           [ 1.36969701e-01, -5.78849752e-01, -7.49388510e-01,
            -8.56650911e-01],
           [ 3.90687328e-02, -1.43062213e+00, -7.18055923e-01,
            -7.71225050e-01],
           [-3.08014162e+00, -3.18109245e-01, -8.74006722e-01,
            -1.43427393e+00],
           [-7.39711134e-01, -9.08422959e-01, -9.67453062e-01,
            -6.72212679e-01],
           [-1.44072036e+00, -2.02622075e+00, -1.34099551e+00,
            -7.04179818e-01],
           [-3.80528374e+00, -5.75288486e-01, -2.76442766e+00,
            -5.10925638e+00],
           [-3.20494027e-02, -1.75219858e+00, -1.37288887e+00,
            -1.17744995e+00]])

  • Yhat
    (observations, response_vars)
    float64
    2.362e+03 1.194e+04 ... 5.448e+05

    array([[  2362.29402198,  11935.07585878,   1496.45841003,
            443789.39246347],
           [  1067.93948791,   7172.26815602,   1884.33993753,
            537793.90966151],
           [   925.44151256,   6975.88715516,   1736.27518851,
            485757.00623504],
           [  1258.38976695,   8152.5027584 ,   1936.75790912,
            525832.89616824],
           [  2282.48175499,  11172.66874738,   1384.55898832,
            397683.04245927],
           [  1129.81865706,   6938.89558614,   1857.37859987,
            530881.21546036],
           [  1003.41968191,   7511.9407956 ,   1949.73685979,
            538564.48326651],
           [   697.51586092,   7384.42010339,   1650.21749434,
            493389.51034526],
           [  1056.24965701,   6008.95776347,   1641.74170094,
            485024.68695605],
           [  1310.64711747,   8028.79601801,   1616.02611257,
            437414.69671848],
    ...
           [  1419.16213464,   6337.30291911,   1679.71467477,
            497225.21001823],
           [  1274.06596048,   6080.04041044,   1622.74714284,
            496286.69338804],
           [  1351.13448897,   5518.25084463,   1461.72297133,
            449733.5352067 ],
           [   686.60743877,   7197.77585286,   1649.37440887,
            497081.83133045],
           [  1131.947039  ,   6679.50205077,   1856.15657721,
            535445.78385057],
           [  1056.53268552,   6381.94697017,   1642.96643645,
            478447.90280334],
           [  1067.27854769,   7115.06206355,   1884.31004692,
            538794.04273634],
           [  3217.37012436,  14612.7997061 ,   1597.00462335,
            476413.78761054],
           [  1394.36749824,   6849.45914991,   1905.98428271,
            500053.5820539 ],
           [  1040.07121356,   7097.70238219,   1922.01751882,
            544844.83934527]])

  • statistics
    (response_vars, statistic)
    float64
    0.3927 0.01852 ... 0.5869 0.995

    array([[ 3.92662580e-01,  1.85185185e-02,  3.47738755e-01,
            -7.07614281e+00,  7.48407335e-01,  3.83796280e-01,
             4.75098802e+02,  5.12638064e-01,  7.12385359e-16,
             6.16203720e-01,  9.80839542e-01],
           [ 2.26325150e-01,  1.51851852e-02,  4.31630931e-01,
            -8.52833615e+00,  1.16216958e+00,  2.26195300e-01,
             3.44042165e+03,  2.62137505e-01,  9.67265863e-05,
             7.73804700e-01,  9.84698136e-01],
           [ 3.02031718e-01,  1.79629630e-02,  9.47306916e-02,
            -5.66326262e+00,  1.21926342e+00,  3.00212585e-01,
             1.97373557e+02,  5.24081050e-01,  1.22903113e-16,
             6.99787415e-01,  9.91782135e-01],
           [ 4.21981577e-01,  2.68518519e-02,  6.34206739e-02,
            -1.11127234e+01,  1.10483669e+00,  4.13094322e-01,
             3.75029928e+04,  6.41002735e-01,  2.17759263e-26,
             5.86905678e-01,  9.95037008e-01]])

• Indexes: (6)
  • observations
    PandasIndex

    PandasIndex(Index([ 756,  769,  692,  616,   35,  164,  680,  331,  299,  727,
           ...
             27,  959,   29,  346,  304,  264,  798,  751,  470, 1043],
          dtype='int64', name='observations', length=216))

  • response_vars
    PandasIndex

    PandasIndex(Index(['WM-hypointensities', 'Right-Lateral-Ventricle', 'Right-Amygdala',
           'CortexVol'],
          dtype='object', name='response_vars'))

  • covariates
    PandasIndex

    PandasIndex(Index(['age'], dtype='object', name='covariates'))

  • batch_effect_dims
    PandasIndex

    PandasIndex(Index(['sex', 'site'], dtype='object', name='batch_effect_dims'))

  • centile
    PandasIndex

    PandasIndex(Index([0.05, 0.25, 0.5, 0.75, 0.95], dtype='float64', name='centile'))

  • statistic
    PandasIndex

    PandasIndex(Index(['EXPV', 'MACE', 'MAPE', 'MSLL', 'NLL', 'R2', 'RMSE', 'Rho', 'Rho_p',
           'SMSE', 'ShapiroW'],
          dtype='object', name='statistic'))

• Attributes: (7)

  real_ids :

  True

  is_scaled :

  False

  name :

  fcon1000_test

  unique_batch_effects :

  {np.str_('sex'): [np.str_('F'), np.str_('M')], np.str_('site'): [np.str_('AnnArbor_a'), np.str_('AnnArbor_b'), np.str_('Atlanta'), np.str_('Baltimore'), np.str_('Bangor'), np.str_('Beijing_Zang'), np.str_('Berlin_Margulies'), np.str_('Cambridge_Buckner'), np.str_('Cleveland'), np.str_('ICBM'), np.str_('Leiden_2180'), np.str_('Leiden_2200'), np.str_('Milwaukee_b'), np.str_('Munchen'), np.str_('NewYork_a'), np.str_('NewYork_a_ADHD'), np.str_('Newark'), np.str_('Oulu'), np.str_('Oxford'), np.str_('PaloAlto'), np.str_('Pittsburgh'), np.str_('Queensland'), np.str_('SaintLouis')]}

  batch_effect_counts :

  defaultdict(<function NormData.register_batch_effects.<locals>.<lambda> at 0x15fefdee0>, {np.str_('sex'): {np.str_('F'): 589, np.str_('M'): 489}, np.str_('site'): {np.str_('AnnArbor_a'): 24, np.str_('AnnArbor_b'): 32, np.str_('Atlanta'): 28, np.str_('Baltimore'): 23, np.str_('Bangor'): 20, np.str_('Beijing_Zang'): 198, np.str_('Berlin_Margulies'): 26, np.str_('Cambridge_Buckner'): 198, np.str_('Cleveland'): 31, np.str_('ICBM'): 85, np.str_('Leiden_2180'): 12, np.str_('Leiden_2200'): 19, np.str_('Milwaukee_b'): 46, np.str_('Munchen'): 15, np.str_('NewYork_a'): 83, np.str_('NewYork_a_ADHD'): 25, np.str_('Newark'): 19, np.str_('Oulu'): 102, np.str_('Oxford'): 22, np.str_('PaloAlto'): 17, np.str_('Pittsburgh'): 3, np.str_('Queensland'): 19, np.str_('SaintLouis'): 31}})

  covariate_ranges :

  {np.str_('age'): {'mean': np.float64(28.251224489795916), 'min': np.float64(7.88), 'max': np.float64(85.0)}}

  batch_effect_covariate_ranges :

  {np.str_('sex'): {np.str_('F'): {np.str_('age'): {'mean': np.float64(28.05332767402377), 'min': np.float64(7.88), 'max': np.float64(85.0)}}, np.str_('M'): {np.str_('age'): {'mean': np.float64(28.48959100204499), 'min': np.float64(9.21), 'max': np.float64(78.0)}}}, np.str_('site'): {np.str_('AnnArbor_a'): {np.str_('age'): {'mean': np.float64(21.28333333333333), 'min': np.float64(13.41), 'max': np.float64(40.98)}}, np.str_('AnnArbor_b'): {np.str_('age'): {'mean': np.float64(44.40625), 'min': np.float64(19.0), 'max': np.float64(79.0)}}, np.str_('Atlanta'): {np.str_('age'): {'mean': np.float64(30.892857142857142), 'min': np.float64(22.0), 'max': np.float64(57.0)}}, np.str_('Baltimore'): {np.str_('age'): {'mean': np.float64(29.26086956521739), 'min': np.float64(20.0), 'max': np.float64(40.0)}}, np.str_('Bangor'): {np.str_('age'): {'mean': np.float64(23.4), 'min': np.float64(19.0), 'max': np.float64(38.0)}}, np.str_('Beijing_Zang'): {np.str_('age'): {'mean': np.float64(21.161616161616163), 'min': np.float64(18.0), 'max': np.float64(26.0)}}, np.str_('Berlin_Margulies'): {np.str_('age'): {'mean': np.float64(29.76923076923077), 'min': np.float64(23.0), 'max': np.float64(44.0)}}, np.str_('Cambridge_Buckner'): {np.str_('age'): {'mean': np.float64(21.03030303030303), 'min': np.float64(18.0), 'max': np.float64(30.0)}}, np.str_('Cleveland'): {np.str_('age'): {'mean': np.float64(43.54838709677419), 'min': np.float64(24.0), 'max': np.float64(60.0)}}, np.str_('ICBM'): {np.str_('age'): {'mean': np.float64(44.04705882352941), 'min': np.float64(19.0), 'max': np.float64(85.0)}}, np.str_('Leiden_2180'): {np.str_('age'): {'mean': np.float64(23.0), 'min': np.float64(20.0), 'max': np.float64(27.0)}}, np.str_('Leiden_2200'): {np.str_('age'): {'mean': np.float64(21.68421052631579), 'min': np.float64(18.0), 'max': np.float64(28.0)}}, np.str_('Milwaukee_b'): {np.str_('age'): {'mean': np.float64(53.58695652173913), 'min': np.float64(44.0), 'max': np.float64(65.0)}}, np.str_('Munchen'): {np.str_('age'): {'mean': np.float64(68.13333333333334), 'min': np.float64(63.0), 'max': np.float64(74.0)}}, np.str_('NewYork_a'): {np.str_('age'): {'mean': np.float64(24.507710843373495), 'min': np.float64(7.88), 'max': np.float64(49.16)}}, np.str_('NewYork_a_ADHD'): {np.str_('age'): {'mean': np.float64(34.9952), 'min': np.float64(20.69), 'max': np.float64(50.9)}}, np.str_('Newark'): {np.str_('age'): {'mean': np.float64(24.105263157894736), 'min': np.float64(21.0), 'max': np.float64(39.0)}}, np.str_('Oulu'): {np.str_('age'): {'mean': np.float64(21.519607843137255), 'min': np.float64(20.0), 'max': np.float64(23.0)}}, np.str_('Oxford'): {np.str_('age'): {'mean': np.float64(29.0), 'min': np.float64(20.0), 'max': np.float64(35.0)}}, np.str_('PaloAlto'): {np.str_('age'): {'mean': np.float64(32.470588235294116), 'min': np.float64(22.0), 'max': np.float64(46.0)}}, np.str_('Pittsburgh'): {np.str_('age'): {'mean': np.float64(32.333333333333336), 'min': np.float64(25.0), 'max': np.float64(47.0)}}, np.str_('Queensland'): {np.str_('age'): {'mean': np.float64(25.94736842105263), 'min': np.float64(20.0), 'max': np.float64(34.0)}}, np.str_('SaintLouis'): {np.str_('age'): {'mean': np.float64(25.096774193548388), 'min': np.float64(21.0), 'max': np.float64(29.0)}}}}

Plot the results

The PCNtoolkit offers are a number of different plotting functions: 1. plot_centiles: Plot the predicted centiles for a model on top of harmonized scatter data. 2. plot_centiles_advanced: a more advanced version of plot_centiles, with more configuration options, coloring, and conditionals 3. plot_qq: Plot the QQ-plot of the predicted Z-scores 4. plot_ridge: Plot density plots of the predicted Z-scores

Let’s start with the centiles.

plot_centiles_advanced(
    model,
    centiles=[0.05, 0.5, 0.95],  # Plot these centiles, the default is [0.05, 0.25, 0.5, 0.75, 0.95]
    scatter_data=train,  # Scatter this data along with the centiles
    batch_effects={"site": ["Beijing_Zang", "AnnArbor_a"], "sex": ["M"]},  # Highlight these groups
    show_other_data=True,  # scatter data not in those groups as smaller black circles
    harmonize_data=True,  # harmonize the scatterdata, this means that we 'remove' the batch effects from the data, by simulating what the data would have looked like if all data was from the same batch.
    conditionals=[30]
)

Now let’s see the qq plots

plot_qq(test, plot_id_line=True)

We can also split the QQ plots by batch effects:

plot_qq(test, plot_id_line=True, hue_data="sex", split_data="sex")
sns.set_theme(style="darkgrid", rc={"axes.facecolor": (0, 0, 0, 0)})

And finally the ridge plot:

plot_ridge(
    train, "Z", split_by="sex"
)  # We can also show the 'Y' variable, and that will show the marginal distribution of the response variable, per batch effect.

/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/pcntoolkit/util/plotter.py:817: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  plt.tight_layout()

/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/pcntoolkit/util/plotter.py:817: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  plt.tight_layout()

/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/pcntoolkit/util/plotter.py:817: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  plt.tight_layout()

/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/seaborn/axisgrid.py:123: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  self._figure.tight_layout(*args, **kwargs)
/opt/anaconda3/envs/ptk/lib/python3.12/site-packages/pcntoolkit/util/plotter.py:817: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.
  plt.tight_layout()

Evaluation statistcs are stored in the NormData object:

display(train.get_statistics_df())
display(test.get_statistics_df())

statistic	EXPV	MACE	MAPE	MSLL	NLL	R2	RMSE	Rho	Rho_p	SMSE	ShapiroW
response_vars
CortexVol	0.519196	0.010719	0.057414	-11.208387	1.072185	0.519179	36152.571241	0.710485	2.199838e-133	0.480821	0.995229
Right-Amygdala	0.383195	0.013689	0.088632	-5.738379	1.181817	0.383064	192.435690	0.599840	2.338456e-85	0.616936	0.995682
Right-Lateral-Ventricle	0.214931	0.013968	0.400007	-8.643762	1.031061	0.213011	3415.620496	0.371987	1.116620e-29	0.786989	0.987490
WM-hypointensities	0.350038	0.014664	0.320848	-7.437633	0.686587	0.349278	658.816151	0.488287	7.601376e-53	0.650722	0.986337

statistic	EXPV	MACE	MAPE	MSLL	NLL	R2	RMSE	Rho	Rho_p	SMSE	ShapiroW
response_vars
CortexVol	0.421982	0.026852	0.063421	-11.112723	1.104837	0.413094	37502.992757	0.641003	2.177593e-26	0.586906	0.995037
Right-Amygdala	0.302032	0.017963	0.094731	-5.663263	1.219263	0.300213	197.373557	0.524081	1.229031e-16	0.699787	0.991782
Right-Lateral-Ventricle	0.226325	0.015185	0.431631	-8.528336	1.162170	0.226195	3440.421655	0.262138	9.672659e-05	0.773805	0.984698
WM-hypointensities	0.392663	0.018519	0.347739	-7.076143	0.748407	0.383796	475.098802	0.512638	7.123854e-16	0.616204	0.980840

What’s next?

Now we have a normative hierarchical Bayesian regression model, we can use it to:

• Make predictions on new data

• Harmonize data, this means that we ‘remove’ the batch effects from the data, by simulating what the data would have looked like if all data was from the same batch.

• Synthesize new data

• Extend the model using data from new batches

Predicting

model.predict(test)

<xarray.NormData> Size: 98kB
Dimensions:            (observations: 216, response_vars: 4, covariates: 1,
                        batch_effect_dims: 2, statistic: 11, centile: 5)
Coordinates:
  * observations       (observations) int64 2kB 756 769 692 616 ... 751 470 1043
  * response_vars      (response_vars) <U23 368B 'WM-hypointensities' ... 'Co...
  * covariates         (covariates) <U3 12B 'age'
  * batch_effect_dims  (batch_effect_dims) <U4 32B 'sex' 'site'
  * statistic          (statistic) <U8 352B 'EXPV' 'MACE' ... 'SMSE' 'ShapiroW'
  * centile            (centile) float64 40B 0.05 0.25 0.5 0.75 0.95
Data variables:
    subject_ids        (observations) object 2kB 'Munchen_sub96752' ... 'Quee...
    Y                  (observations, response_vars) float64 7kB 2.721e+03 .....
    X                  (observations, covariates) float64 2kB 63.0 ... 23.0
    batch_effects      (observations, batch_effect_dims) <U17 29kB 'F' ... 'Q...
    Z                  (observations, response_vars) float64 7kB 0.5221 ... -...
    logp               (observations, response_vars) float64 7kB -1.623 ... -...
    Yhat               (observations, response_vars) float64 7kB 2.362e+03 .....
    statistics         (response_vars, statistic) float64 352B 0.3927 ... 0.995
    centiles           (centile, observations, response_vars) float64 35kB 1....
Attributes:
    real_ids:                       True
    is_scaled:                      False
    name:                           fcon1000_test
    unique_batch_effects:           {np.str_('sex'): [np.str_('F'), np.str_('...
    batch_effect_counts:            defaultdict(<function NormData.register_b...
    covariate_ranges:               {np.str_('age'): {'mean': np.float64(28.2...
    batch_effect_covariate_ranges:  {np.str_('sex'): {np.str_('F'): {np.str_(...

xarray.NormData

• Dimensions:
  • observations: 216
  • response_vars: 4
  • covariates: 1
  • batch_effect_dims: 2
  • statistic: 11
  • centile: 5
• Coordinates: (6)
  • observations
    (observations)
    int64
    756 769 692 616 ... 751 470 1043

    array([ 756,  769,  692, ...,  751,  470, 1043], shape=(216,))

  • response_vars
    (response_vars)
    <U23
    'WM-hypointensities' ... 'Cortex...

    array(['WM-hypointensities', 'Right-Lateral-Ventricle', 'Right-Amygdala',
           'CortexVol'], dtype='<U23')

  • covariates
    (covariates)
    <U3
    'age'

    array(['age'], dtype='<U3')

  • batch_effect_dims
    (batch_effect_dims)
    <U4
    'sex' 'site'

    array(['sex', 'site'], dtype='<U4')

  • statistic
    (statistic)
    <U8
    'EXPV' 'MACE' ... 'SMSE' 'ShapiroW'

    array(['EXPV', 'MACE', 'MAPE', 'MSLL', 'NLL', 'R2', 'RMSE', 'Rho', 'Rho_p',
           'SMSE', 'ShapiroW'], dtype='<U8')

  • centile
    (centile)
    float64
    0.05 0.25 0.5 0.75 0.95

    array([0.05, 0.25, 0.5 , 0.75, 0.95])

• Data variables: (9)
  • subject_ids
    (observations)
    object
    'Munchen_sub96752' ... 'Queensla...

    array(['Munchen_sub96752', 'NewYork_a_sub18638', 'Leiden_2200_sub87320',
           'ICBM_sub47658', 'AnnArbor_b_sub45569', 'Beijing_Zang_sub18960',
           'Leiden_2200_sub18456', 'Berlin_Margulies_sub27711',
           'Beijing_Zang_sub87776', 'Milwaukee_b_sub63196',
           'Beijing_Zang_sub07144', 'Atlanta_sub76280',
           'Beijing_Zang_sub40037', 'Cambridge_Buckner_sub17737',
           'ICBM_sub89049', 'ICBM_sub55656', 'Oulu_sub45566',
           'Beijing_Zang_sub89088', 'Atlanta_sub16563',
           'Cambridge_Buckner_sub51172', 'Oulu_sub98739',
           'Queensland_sub49845', 'Cambridge_Buckner_sub84256',
           'Cleveland_sub80263', 'ICBM_sub16607', 'Newark_sub46570',
           'NewYork_a_sub88286', 'Cambridge_Buckner_sub02591',
           'Oulu_sub66467', 'Beijing_Zang_sub74386', 'Newark_sub55760',
           'ICBM_sub30623', 'Oulu_sub68752', 'Leiden_2180_sub19281',
           'Beijing_Zang_sub50972', 'Beijing_Zang_sub85030',
           'Milwaukee_b_sub36386', 'Baltimore_sub31837', 'PaloAlto_sub84978',
           'Oulu_sub01077', 'NewYork_a_ADHD_sub54828', 'PaloAlto_sub96705',
           'Cambridge_Buckner_sub40635', 'ICBM_sub66794',
           'Beijing_Zang_sub46541', 'Beijing_Zang_sub87089',
           'Pittsburgh_sub97823', 'Beijing_Zang_sub98617', 'ICBM_sub92028',
    ...
           'Leiden_2200_sub04484', 'Beijing_Zang_sub80163', 'ICBM_sub02382',
           'Cambridge_Buckner_sub77435', 'NewYork_a_sub54887',
           'Oulu_sub85532', 'Baltimore_sub73823', 'Beijing_Zang_sub29590',
           'Oulu_sub99718', 'Beijing_Zang_sub08455', 'Beijing_Zang_sub85543',
           'Cambridge_Buckner_sub45354', 'Beijing_Zang_sub07717',
           'Baltimore_sub76160', 'Beijing_Zang_sub17093',
           'AnnArbor_b_sub90127', 'SaintLouis_sub73002',
           'Queensland_sub93238', 'Cleveland_sub34189',
           'Cambridge_Buckner_sub89107', 'Atlanta_sub75153',
           'NewYork_a_ADHD_sub73035', 'Cambridge_Buckner_sub59434',
           'Milwaukee_b_sub44912', 'Cleveland_sub46739', 'Oulu_sub20495',
           'SaintLouis_sub28304', 'Cambridge_Buckner_sub35430',
           'Oulu_sub86362', 'Newark_sub58526', 'Leiden_2180_sub12255',
           'ICBM_sub48210', 'Cambridge_Buckner_sub77989',
           'Berlin_Margulies_sub75506', 'NewYork_a_sub29216',
           'Beijing_Zang_sub05267', 'AnnArbor_b_sub18546', 'Oulu_sub75620',
           'AnnArbor_b_sub30250', 'Berlin_Margulies_sub86111',
           'Beijing_Zang_sub89592', 'Beijing_Zang_sub68012',
           'NewYork_a_sub50559', 'Munchen_sub66933',
           'Cambridge_Buckner_sub59729', 'Queensland_sub86245'], dtype=object)

  • Y
    (observations, response_vars)
    float64
    2.721e+03 1.289e+04 ... 5.035e+05

    array([[2.72140000e+03, 1.28916000e+04, 1.43940000e+03, 4.57858328e+05],
           [1.14310000e+03, 9.91910000e+03, 1.64970000e+03, 5.26780362e+05],
           [9.55800000e+02, 7.47730000e+03, 1.83850000e+03, 4.95744471e+05],
           [1.47390000e+03, 1.43021000e+04, 1.86770000e+03, 5.85303839e+05],
           [7.57800000e+02, 4.11930000e+03, 1.32500000e+03, 3.33111552e+05],
           [8.71100000e+02, 5.03090000e+03, 1.90730000e+03, 5.10794940e+05],
           [1.20730000e+03, 1.78664000e+04, 2.02220000e+03, 5.50533325e+05],
           [5.95000000e+02, 5.00790000e+03, 2.01070000e+03, 4.67673977e+05],
           [6.82400000e+02, 7.28660000e+03, 1.45630000e+03, 4.60129533e+05],
           [4.45100000e+02, 5.74290000e+03, 1.47450000e+03, 4.44494817e+05],
           [1.62000000e+03, 3.71370000e+03, 2.00110000e+03, 5.59424624e+05],
           [6.02800000e+02, 5.30120000e+03, 1.36100000e+03, 4.21551234e+05],
           [1.43250000e+03, 4.42970000e+03, 1.65080000e+03, 5.19842763e+05],
           [1.90820000e+03, 3.57810000e+03, 1.88370000e+03, 5.06679262e+05],
           [1.83400000e+03, 3.27190000e+03, 2.05120000e+03, 5.35569987e+05],
           [4.59600000e+02, 3.98580000e+03, 1.45470000e+03, 4.67607555e+05],
           [1.21000000e+03, 8.72130000e+03, 1.71430000e+03, 5.30904612e+05],
           [8.45900000e+02, 6.59310000e+03, 1.61900000e+03, 5.09371867e+05],
           [9.95200000e+02, 7.04020000e+03, 1.99490000e+03, 4.60068379e+05],
           [1.73470000e+03, 4.01480000e+03, 1.51620000e+03, 4.87269373e+05],
    ...
           [7.85800000e+02, 5.70900000e+03, 1.47480000e+03, 4.53982166e+05],
           [2.24010000e+03, 4.36660000e+03, 2.04210000e+03, 5.58453123e+05],
           [7.58100000e+02, 6.52980000e+03, 1.56730000e+03, 4.73575183e+05],
           [1.44050000e+03, 6.70530000e+03, 1.20540000e+03, 3.82788491e+05],
           [8.18600000e+02, 9.38330000e+03, 1.96740000e+03, 5.02713911e+05],
           [3.76990000e+03, 1.58644000e+04, 1.79170000e+03, 5.12490348e+05],
           [8.80200000e+02, 4.37020000e+03, 1.75520000e+03, 4.37300069e+05],
           [8.23900000e+02, 6.37900000e+03, 1.57650000e+03, 5.67331908e+05],
           [2.11390000e+03, 1.07225000e+04, 1.84380000e+03, 5.12273764e+05],
           [7.41900000e+02, 8.80170000e+03, 1.60640000e+03, 4.91973562e+05],
           [1.33390000e+03, 6.98000000e+03, 1.74850000e+03, 4.78907154e+05],
           [7.07300000e+02, 5.68070000e+03, 1.53450000e+03, 4.74077083e+05],
           [1.13410000e+03, 5.59220000e+03, 1.62620000e+03, 4.54163909e+05],
           [4.38600000e+02, 6.33000000e+03, 1.59670000e+03, 4.68067037e+05],
           [9.66300000e+02, 9.21550000e+03, 1.78250000e+03, 5.09199708e+05],
           [4.24300000e+02, 4.51110000e+03, 1.70200000e+03, 5.26635258e+05],
           [6.04700000e+02, 7.59080000e+03, 1.69930000e+03, 5.20499663e+05],
           [2.34320000e+03, 1.71923000e+04, 1.79380000e+03, 4.86680791e+05],
           [2.72170000e+03, 6.08600000e+03, 2.32470000e+03, 6.10402006e+05],
           [7.03500000e+02, 1.07003000e+04, 1.67620000e+03, 5.03535771e+05]])

  • X
    (observations, covariates)
    float64
    63.0 23.27 22.0 ... 72.0 23.0 23.0

    array([[63.  ],
           [23.27],
           [22.  ],
           [42.  ],
           [63.  ],
           [23.  ],
           [21.  ],
           [26.  ],
           [21.  ],
           [49.  ],
           [20.  ],
           [23.  ],
           [20.  ],
           [26.  ],
           [35.  ],
           [21.  ],
           [22.  ],
           [19.  ],
           [34.  ],
           [18.  ],
    ...
           [21.  ],
           [20.  ],
           [22.  ],
           [25.  ],
           [25.  ],
           [73.  ],
           [22.  ],
           [28.  ],
           [29.06],
           [19.  ],
           [20.  ],
           [22.  ],
           [19.  ],
           [24.  ],
           [21.  ],
           [24.  ],
           [22.79],
           [72.  ],
           [23.  ],
           [23.  ]])

  • batch_effects
    (observations, batch_effect_dims)
    <U17
    'F' 'Munchen' ... 'M' 'Queensland'

    array([['F', 'Munchen'],
           ['M', 'NewYork_a'],
           ['F', 'Leiden_2200'],
           ['M', 'ICBM'],
           ['F', 'AnnArbor_b'],
           ['M', 'Beijing_Zang'],
           ['M', 'Leiden_2200'],
           ['F', 'Berlin_Margulies'],
           ['F', 'Beijing_Zang'],
           ['F', 'Milwaukee_b'],
           ['M', 'Beijing_Zang'],
           ['F', 'Atlanta'],
           ['F', 'Beijing_Zang'],
           ['F', 'Cambridge_Buckner'],
           ['M', 'ICBM'],
           ['F', 'ICBM'],
           ['M', 'Oulu'],
           ['F', 'Beijing_Zang'],
           ['M', 'Atlanta'],
           ['F', 'Cambridge_Buckner'],
    ...
           ['F', 'SaintLouis'],
           ['M', 'Cambridge_Buckner'],
           ['F', 'Oulu'],
           ['F', 'Newark'],
           ['M', 'Leiden_2180'],
           ['M', 'ICBM'],
           ['F', 'Cambridge_Buckner'],
           ['M', 'Berlin_Margulies'],
           ['M', 'NewYork_a'],
           ['F', 'Beijing_Zang'],
           ['M', 'AnnArbor_b'],
           ['F', 'Oulu'],
           ['F', 'AnnArbor_b'],
           ['F', 'Berlin_Margulies'],
           ['M', 'Beijing_Zang'],
           ['F', 'Beijing_Zang'],
           ['M', 'NewYork_a'],
           ['M', 'Munchen'],
           ['M', 'Cambridge_Buckner'],
           ['M', 'Queensland']], dtype='<U17')

  • Z
    (observations, response_vars)
    float64
    0.5221 0.3941 ... -1.314 -1.141

    array([[ 5.22111034e-01,  3.94099738e-01, -2.66792814e-01,
             4.06460635e-01],
           [ 4.22466208e-01,  9.65448035e-01, -1.23115188e+00,
            -2.75083913e-01],
           [ 3.18291793e-01,  4.04732097e-01,  5.45620314e-01,
             2.88865481e-01],
           [ 5.69111007e-01,  1.36465454e+00, -2.58229393e-01,
             1.48401663e+00],
           [-2.02047143e+00, -1.70824599e+00, -2.79573577e-01,
            -1.90516036e+00],
           [-5.72758421e-01, -5.16773934e-01,  3.06331316e-01,
            -5.25596222e-01],
           [ 6.75159082e-01,  2.69725411e+00,  4.09408581e-01,
             3.40759379e-01],
           [-2.58377585e-02, -7.19266423e-01,  1.68705934e+00,
            -6.77610473e-01],
           [-1.09910335e+00,  6.32661043e-01, -9.32247384e-01,
            -6.61896768e-01],
           [-1.64782803e+00, -3.09579804e-01, -6.85539790e-01,
             2.02469096e-01],
    ...
           [ 3.23383225e-02,  4.56707221e-01,  3.94488157e-01,
            -4.75028151e-01],
           [-2.28010369e+00,  1.25793914e-01, -3.72871224e-01,
            -5.86864722e-01],
           [-3.68727350e-01,  2.77777856e-01,  8.28827400e-01,
             1.39496191e-01],
           [-5.46991665e-01, -3.98633468e-02, -1.89215774e-01,
            -7.79284242e-01],
           [-2.02910739e-01,  9.55953072e-01, -2.96613602e-01,
            -7.00202583e-01],
           [-2.75327977e+00, -4.66405172e-01,  3.47307354e-01,
             1.29650564e+00],
           [-1.66030307e+00,  3.97481554e-01, -9.24188289e-01,
            -4.76654211e-01],
           [-3.08981685e-01,  5.44325984e-01,  9.89606313e-01,
             2.91387578e-01],
           [ 2.46018649e+00,  2.29349516e-03,  1.98860470e+00,
             3.02005316e+00],
           [-9.65105272e-01,  1.16469286e+00, -1.31411301e+00,
            -1.14085260e+00]])

  • logp
    (observations, response_vars)
    float64
    -1.623 -1.648 ... -1.373 -1.177

    array([[-1.62282548e+00, -1.64785626e+00, -6.29004113e-01,
            -7.04684293e-01],
           [-4.37893961e-01, -1.50704043e+00, -1.24324841e+00,
            -6.16205101e-01],
           [-3.53409858e-01, -8.94472096e-01, -9.82346850e-01,
            -6.67329860e-01],
           [-9.91335914e-01, -2.35266878e+00, -7.68952127e-01,
            -1.77734024e+00],
           [-2.46130899e+00, -1.80307654e+00, -6.15292810e-01,
            -2.29610528e+00],
           [ 1.69843795e-01, -2.76843190e-01, -8.49856792e-01,
            -6.82519637e-01],
           [-7.19451109e-01, -4.77541627e+00, -9.05046781e-01,
            -6.85689292e-01],
           [-6.04283284e-02, -3.85498427e-01, -2.26037994e+00,
            -8.00264151e-01],
           [-5.70768126e-02, -1.08106993e+00, -9.50155776e-01,
            -7.47096740e-01],
           [-1.34655206e+00, -7.91392769e-01, -8.05409464e-01,
            -6.93195709e-01],
    ...
           [-1.44179101e-01, -8.85692449e-01, -8.92070483e-01,
            -6.88624281e-01],
           [-1.90784482e+00, -6.47669808e-01, -7.26611956e-01,
            -7.10456058e-01],
           [ 4.77865405e-02, -6.97197063e-01, -1.18036135e+00,
            -6.40334051e-01],
           [ 1.36969701e-01, -5.78849752e-01, -7.49388510e-01,
            -8.56650911e-01],
           [ 3.90687328e-02, -1.43062213e+00, -7.18055923e-01,
            -7.71225050e-01],
           [-3.08014162e+00, -3.18109245e-01, -8.74006722e-01,
            -1.43427393e+00],
           [-7.39711134e-01, -9.08422959e-01, -9.67453062e-01,
            -6.72212679e-01],
           [-1.44072036e+00, -2.02622075e+00, -1.34099551e+00,
            -7.04179818e-01],
           [-3.80528374e+00, -5.75288486e-01, -2.76442766e+00,
            -5.10925638e+00],
           [-3.20494027e-02, -1.75219858e+00, -1.37288887e+00,
            -1.17744995e+00]])

  • Yhat
    (observations, response_vars)
    float64
    2.362e+03 1.194e+04 ... 5.448e+05

    array([[  2362.29402198,  11935.07585878,   1496.45841003,
            443789.39246347],
           [  1067.93948791,   7172.26815602,   1884.33993753,
            537793.90966151],
           [   925.44151256,   6975.88715516,   1736.27518851,
            485757.00623504],
           [  1258.38976695,   8152.5027584 ,   1936.75790912,
            525832.89616824],
           [  2282.48175499,  11172.66874738,   1384.55898832,
            397683.04245927],
           [  1129.81865706,   6938.89558614,   1857.37859987,
            530881.21546036],
           [  1003.41968191,   7511.9407956 ,   1949.73685979,
            538564.48326651],
           [   697.51586092,   7384.42010339,   1650.21749434,
            493389.51034526],
           [  1056.24965701,   6008.95776347,   1641.74170094,
            485024.68695605],
           [  1310.64711747,   8028.79601801,   1616.02611257,
            437414.69671848],
    ...
           [  1419.16213464,   6337.30291911,   1679.71467477,
            497225.21001823],
           [  1274.06596048,   6080.04041044,   1622.74714284,
            496286.69338804],
           [  1351.13448897,   5518.25084463,   1461.72297133,
            449733.5352067 ],
           [   686.60743877,   7197.77585286,   1649.37440887,
            497081.83133045],
           [  1131.947039  ,   6679.50205077,   1856.15657721,
            535445.78385057],
           [  1056.53268552,   6381.94697017,   1642.96643645,
            478447.90280334],
           [  1067.27854769,   7115.06206355,   1884.31004692,
            538794.04273634],
           [  3217.37012436,  14612.7997061 ,   1597.00462335,
            476413.78761054],
           [  1394.36749824,   6849.45914991,   1905.98428271,
            500053.5820539 ],
           [  1040.07121356,   7097.70238219,   1922.01751882,
            544844.83934527]])

  • statistics
    (response_vars, statistic)
    float64
    0.3927 0.01852 ... 0.5869 0.995

    array([[ 3.92662580e-01,  1.85185185e-02,  3.47738755e-01,
            -7.07614281e+00,  7.48407335e-01,  3.83796280e-01,
             4.75098802e+02,  5.12638064e-01,  7.12385359e-16,
             6.16203720e-01,  9.80839542e-01],
           [ 2.26325150e-01,  1.51851852e-02,  4.31630931e-01,
            -8.52833615e+00,  1.16216958e+00,  2.26195300e-01,
             3.44042165e+03,  2.62137505e-01,  9.67265863e-05,
             7.73804700e-01,  9.84698136e-01],
           [ 3.02031718e-01,  1.79629630e-02,  9.47306916e-02,
            -5.66326262e+00,  1.21926342e+00,  3.00212585e-01,
             1.97373557e+02,  5.24081050e-01,  1.22903113e-16,
             6.99787415e-01,  9.91782135e-01],
           [ 4.21981577e-01,  2.68518519e-02,  6.34206739e-02,
            -1.11127234e+01,  1.10483669e+00,  4.13094322e-01,
             3.75029928e+04,  6.41002735e-01,  2.17759263e-26,
             5.86905678e-01,  9.95037008e-01]])

  • centiles
    (centile, observations, response_vars)
    float64
    1.004e+03 4.959e+03 ... 6.059e+05

    array([[[1.00398750e+03, 4.95918729e+03, 1.23858886e+03, 3.86736055e+05],
            [6.08602095e+02, 3.81136404e+03, 1.58699256e+03, 4.79606192e+05],
            [4.61188041e+02, 3.73382051e+03, 1.44083473e+03, 4.27724793e+05],
            ...,
            [1.09822390e+03, 5.73774119e+03, 1.30038803e+03, 4.16559291e+05],
            [9.34238404e+02, 3.51395053e+03, 1.60910669e+03, 4.41917504e+05],
            [5.79942119e+02, 3.76219376e+03, 1.62513992e+03, 4.86708761e+05]],
    
           [[1.52906096e+03, 7.38744766e+03, 1.37175546e+03, 4.18879836e+05],
            [7.86085048e+02, 4.98137572e+03, 1.74054148e+03, 5.12383721e+05],
            [6.40568118e+02, 4.86249913e+03, 1.59339913e+03, 4.60415329e+05],
            ...,
            [1.91753230e+03, 8.82741884e+03, 1.45358310e+03, 4.50275711e+05],
            [1.11202673e+03, 4.67512909e+03, 1.76241305e+03, 4.74666069e+05],
            [7.57730442e+02, 4.92337233e+03, 1.77844629e+03, 5.19457326e+05]],
    
           [[2.07580988e+03, 1.03198971e+04, 1.48322514e+03, 4.42912966e+05],
            [9.70933244e+02, 6.39440744e+03, 1.86911066e+03, 5.36897209e+05],
            [8.27397246e+02, 6.22553695e+03, 1.72114396e+03, 4.84863169e+05],
            ...,
            [2.77052962e+03, 1.25581568e+04, 1.58180962e+03, 4.75490524e+05],
            [1.29719402e+03, 6.07747726e+03, 1.89077914e+03, 4.99157770e+05],
            [9.42897736e+02, 6.32572049e+03, 1.90681237e+03, 5.43949027e+05]],
    
           [[2.94786482e+03, 1.51168517e+04, 1.61037303e+03, 4.67976306e+05],
            [1.26591979e+03, 8.70543430e+03, 2.01572517e+03, 5.62464780e+05],
            [1.12554392e+03, 8.45480271e+03, 1.86681796e+03, 5.10361722e+05],
            ...,
            [4.13092024e+03, 1.86606901e+04, 1.72803919e+03, 5.01790603e+05],
            [1.59268967e+03, 8.37103215e+03, 2.03716196e+03, 5.24702508e+05],
            [1.23839338e+03, 8.61927538e+03, 2.05319520e+03, 5.69493765e+05]],
    
           [[4.65019278e+03, 2.41226868e+04, 1.79662568e+03, 5.03652340e+05],
            [1.84205824e+03, 1.30430486e+04, 2.23036251e+03, 5.98856188e+05],
            [1.70784407e+03, 1.26390713e+04, 2.08007745e+03, 5.46654607e+05],
            ...,
            [6.78643678e+03, 3.01176185e+04, 1.94218858e+03, 5.39228041e+05],
            [2.16982050e+03, 1.26758769e+04, 2.25145994e+03, 5.61061376e+05],
            [1.81552422e+03, 1.29241201e+04, 2.26749317e+03, 6.05852634e+05]]],
          shape=(5, 216, 4))

• Indexes: (6)
  • observations
    PandasIndex

    PandasIndex(Index([ 756,  769,  692,  616,   35,  164,  680,  331,  299,  727,
           ...
             27,  959,   29,  346,  304,  264,  798,  751,  470, 1043],
          dtype='int64', name='observations', length=216))

  • response_vars
    PandasIndex

    PandasIndex(Index(['WM-hypointensities', 'Right-Lateral-Ventricle', 'Right-Amygdala',
           'CortexVol'],
          dtype='object', name='response_vars'))

  • covariates
    PandasIndex

    PandasIndex(Index(['age'], dtype='object', name='covariates'))

  • batch_effect_dims
    PandasIndex

    PandasIndex(Index(['sex', 'site'], dtype='object', name='batch_effect_dims'))

  • statistic
    PandasIndex

    PandasIndex(Index(['EXPV', 'MACE', 'MAPE', 'MSLL', 'NLL', 'R2', 'RMSE', 'Rho', 'Rho_p',
           'SMSE', 'ShapiroW'],
          dtype='object', name='statistic'))

  • centile
    PandasIndex

    PandasIndex(Index([0.05, 0.25, 0.5, 0.75, 0.95], dtype='float64', name='centile'))

• Attributes: (7)

  real_ids :

  True

  is_scaled :

  False

  name :

  fcon1000_test

  unique_batch_effects :

  {np.str_('sex'): [np.str_('F'), np.str_('M')], np.str_('site'): [np.str_('AnnArbor_a'), np.str_('AnnArbor_b'), np.str_('Atlanta'), np.str_('Baltimore'), np.str_('Bangor'), np.str_('Beijing_Zang'), np.str_('Berlin_Margulies'), np.str_('Cambridge_Buckner'), np.str_('Cleveland'), np.str_('ICBM'), np.str_('Leiden_2180'), np.str_('Leiden_2200'), np.str_('Milwaukee_b'), np.str_('Munchen'), np.str_('NewYork_a'), np.str_('NewYork_a_ADHD'), np.str_('Newark'), np.str_('Oulu'), np.str_('Oxford'), np.str_('PaloAlto'), np.str_('Pittsburgh'), np.str_('Queensland'), np.str_('SaintLouis')]}

  batch_effect_counts :

  defaultdict(<function NormData.register_batch_effects.<locals>.<lambda> at 0x15fefdee0>, {np.str_('sex'): {np.str_('F'): 589, np.str_('M'): 489}, np.str_('site'): {np.str_('AnnArbor_a'): 24, np.str_('AnnArbor_b'): 32, np.str_('Atlanta'): 28, np.str_('Baltimore'): 23, np.str_('Bangor'): 20, np.str_('Beijing_Zang'): 198, np.str_('Berlin_Margulies'): 26, np.str_('Cambridge_Buckner'): 198, np.str_('Cleveland'): 31, np.str_('ICBM'): 85, np.str_('Leiden_2180'): 12, np.str_('Leiden_2200'): 19, np.str_('Milwaukee_b'): 46, np.str_('Munchen'): 15, np.str_('NewYork_a'): 83, np.str_('NewYork_a_ADHD'): 25, np.str_('Newark'): 19, np.str_('Oulu'): 102, np.str_('Oxford'): 22, np.str_('PaloAlto'): 17, np.str_('Pittsburgh'): 3, np.str_('Queensland'): 19, np.str_('SaintLouis'): 31}})

  covariate_ranges :

  {np.str_('age'): {'mean': np.float64(28.251224489795916), 'min': np.float64(7.88), 'max': np.float64(85.0)}}

  batch_effect_covariate_ranges :

  {np.str_('sex'): {np.str_('F'): {np.str_('age'): {'mean': np.float64(28.05332767402377), 'min': np.float64(7.88), 'max': np.float64(85.0)}}, np.str_('M'): {np.str_('age'): {'mean': np.float64(28.48959100204499), 'min': np.float64(9.21), 'max': np.float64(78.0)}}}, np.str_('site'): {np.str_('AnnArbor_a'): {np.str_('age'): {'mean': np.float64(21.28333333333333), 'min': np.float64(13.41), 'max': np.float64(40.98)}}, np.str_('AnnArbor_b'): {np.str_('age'): {'mean': np.float64(44.40625), 'min': np.float64(19.0), 'max': np.float64(79.0)}}, np.str_('Atlanta'): {np.str_('age'): {'mean': np.float64(30.892857142857142), 'min': np.float64(22.0), 'max': np.float64(57.0)}}, np.str_('Baltimore'): {np.str_('age'): {'mean': np.float64(29.26086956521739), 'min': np.float64(20.0), 'max': np.float64(40.0)}}, np.str_('Bangor'): {np.str_('age'): {'mean': np.float64(23.4), 'min': np.float64(19.0), 'max': np.float64(38.0)}}, np.str_('Beijing_Zang'): {np.str_('age'): {'mean': np.float64(21.161616161616163), 'min': np.float64(18.0), 'max': np.float64(26.0)}}, np.str_('Berlin_Margulies'): {np.str_('age'): {'mean': np.float64(29.76923076923077), 'min': np.float64(23.0), 'max': np.float64(44.0)}}, np.str_('Cambridge_Buckner'): {np.str_('age'): {'mean': np.float64(21.03030303030303), 'min': np.float64(18.0), 'max': np.float64(30.0)}}, np.str_('Cleveland'): {np.str_('age'): {'mean': np.float64(43.54838709677419), 'min': np.float64(24.0), 'max': np.float64(60.0)}}, np.str_('ICBM'): {np.str_('age'): {'mean': np.float64(44.04705882352941), 'min': np.float64(19.0), 'max': np.float64(85.0)}}, np.str_('Leiden_2180'): {np.str_('age'): {'mean': np.float64(23.0), 'min': np.float64(20.0), 'max': np.float64(27.0)}}, np.str_('Leiden_2200'): {np.str_('age'): {'mean': np.float64(21.68421052631579), 'min': np.float64(18.0), 'max': np.float64(28.0)}}, np.str_('Milwaukee_b'): {np.str_('age'): {'mean': np.float64(53.58695652173913), 'min': np.float64(44.0), 'max': np.float64(65.0)}}, np.str_('Munchen'): {np.str_('age'): {'mean': np.float64(68.13333333333334), 'min': np.float64(63.0), 'max': np.float64(74.0)}}, np.str_('NewYork_a'): {np.str_('age'): {'mean': np.float64(24.507710843373495), 'min': np.float64(7.88), 'max': np.float64(49.16)}}, np.str_('NewYork_a_ADHD'): {np.str_('age'): {'mean': np.float64(34.9952), 'min': np.float64(20.69), 'max': np.float64(50.9)}}, np.str_('Newark'): {np.str_('age'): {'mean': np.float64(24.105263157894736), 'min': np.float64(21.0), 'max': np.float64(39.0)}}, np.str_('Oulu'): {np.str_('age'): {'mean': np.float64(21.519607843137255), 'min': np.float64(20.0), 'max': np.float64(23.0)}}, np.str_('Oxford'): {np.str_('age'): {'mean': np.float64(29.0), 'min': np.float64(20.0), 'max': np.float64(35.0)}}, np.str_('PaloAlto'): {np.str_('age'): {'mean': np.float64(32.470588235294116), 'min': np.float64(22.0), 'max': np.float64(46.0)}}, np.str_('Pittsburgh'): {np.str_('age'): {'mean': np.float64(32.333333333333336), 'min': np.float64(25.0), 'max': np.float64(47.0)}}, np.str_('Queensland'): {np.str_('age'): {'mean': np.float64(25.94736842105263), 'min': np.float64(20.0), 'max': np.float64(34.0)}}, np.str_('SaintLouis'): {np.str_('age'): {'mean': np.float64(25.096774193548388), 'min': np.float64(21.0), 'max': np.float64(29.0)}}}}

Harmonize

# Harmonizing is also easy:
reference_batch_effect = {
    "site": "Beijing_Zang",
    "sex": "M",
}  # Set a pseudo-batch effect. I.e., this means 'pretend that all data was from this site and sex'

# model.harmonize(train, reference_batch_effect=reference_batch_effect)  # <- easy

plt.style.use("seaborn-v0_8")
df = train.to_dataframe()
fig, ax = plt.subplots(1, 2, figsize=(13, 5), sharey=True)
sns.scatterplot(data=df, x=("X", "age"), y=("Y", feature_to_plot), hue=("batch_effects", "site"), ax=ax[0])
sns.scatterplot(data=df, x=("X", "age"), y=("Y_harmonized", feature_to_plot), hue=("batch_effects", "site"), ax=ax[1])
ax[0].title.set_text("Unharmonized")
ax[1].title.set_text("Harmonized")
ax[0].legend([], [])
ax[1].legend([], [])
ax[0].set_xlabel("Age")
ax[0].set_ylabel(feature_to_plot)
ax[1].set_xlabel("Age")
ax[1].set_ylabel(feature_to_plot)
plt.tight_layout()
plt.show()

Synthesize

Our models can synthesize new data that follows the learned distribution.

Not only the distribution of the response variables given a covariate is learned, but also the ranges of the covariates within each batch effect. So if we have fitted a model on a number of sites, and subjects from A have an age between 10 and 20, then the synthesized pseudo-subjects from site A will also have an age between 10 and 20.

Not only that, but we also sample the batch effects in the frequency of the batch effects in the original data. So if the train data contained twice as many subjects from site A as site B, then the synthesized pseudo-subjects will also have twice as many subjects from site A as site B.

# Generate 10000 synthetic datapoints from scratch
synthetic_data = model.synthesize(covariate_range_per_batch_effect=True, n_samples=1000)  # <- also easy
# Show the synthetic data along with the centiles
plot_centiles_advanced(
    model,
    covariate="age",  # Which covariate to plot on the x-axis
    scatter_data=synthetic_data,
)

# Synthesize new Y data for existing X data
new_test_data = test.copy()

# Remove the Y data, this way we will synthesize new Y data for the existing X data
if hasattr(new_test_data, "Y"):
    del new_test_data["Y"]

synthetic = model.synthesize(new_test_data)  # <- will fill in the missing Y data

plot_centiles_advanced(
    model,
    centiles=[0.05, 0.5, 0.95],  # Plot arbitrary centiles
    covariate="age",  # Which covariate to plot on the x-axis
    scatter_data=synthetic,  # Scatter the train data points
    batch_effects="all",  # You can set this to "all" to show all batch effects
    show_other_data=False,  # Show data points that do not match any batch effects
    show_centile_labels=True,
    harmonize_data=True,  # Set this to False to see the difference
    show_legend=False,  # Don't show the legend because it crowds the plot
)

Next steps

Please see the other tutorials for more examples, and we also recommend you to read the documentation! As this toolkit is still in development, the documentation may not be up to date. If you find any issues, please let us know!

Also, feel free to contact us on Github if you have any questions or suggestions.

Have fun modeling!

Bonus content

Here is another model configuration using a SHASH likelihood, but this one also has a linear regression in epsilon and delta. If you have a feature that is heavily skewed and for which the skewness also changes with the covariates, this is the model for you:

# Here's a model with a SHASHb likelihood, with a linear regression in all four parameters, so including epsilon and delta.
# This is a very flexible model, but it will also take a lot longer to run.
mu = make_prior(
    linear=True,
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 10.0)),
    intercept=make_prior(
        random=True,
        mu=make_prior(dist_name="Normal", dist_params=(1.0, 1.0)),
        sigma=make_prior(dist_name="Gamma", dist_params=(3.0, 1.0)),
    ),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
)
sigma = make_prior(
    linear=True,
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 2.0)),
    intercept=make_prior(dist_name="Normal", dist_params=(1.0, 1.0)),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
    mapping="softplus",
    mapping_params=(0.0, 3.0),
)

epsilon = make_prior(
    linear=True,
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 1.0)),
    intercept=make_prior(dist_name="Normal", dist_params=(0.0, 1.0)),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
)

delta = make_prior(
    linear=True,
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 1.0)),
    intercept=make_prior(dist_name="Normal", dist_params=(1.0, 1.0)),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
    mapping="softplus",
    mapping_params=(
        0.0,
        3.0,  # Scale for smoothness
        0.6,  # We need to provide a vertical shift as well, because the SHASH mapping goes a bit wild with low values for delta
    ),
)

shashb2_regression_model = HBR(
    name="template",
    cores=16,
    progressbar=True,
    draws=1500,
    tune=500,
    chains=4,
    nuts_sampler="nutpie",
    likelihood=SHASHbLikelihood(mu, sigma, epsilon, delta),
)