HBR with Beta 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 Beta 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,
    BetaLikelihood,
    make_prior,
    plot_centiles_advanced,
    plot_qq,
    plot_ridge,
)

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

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 BetaLikelihood, which we will use to model our response variable.

# The beta likelihood is a bit different from the other likelihoods, it assumes the response variable is bounded between 0 and 1. The means that we have to apply a min max scaling to the response variable before fitting the model. Luckily, this is built into the PCNtoolkit, we just have to set 'outscaler' to 'minmax' in the `NormativeModel` constructor.

# The beta likelihood takes two parameters, alpha and beta. Intuitively, these encode 'pseudo-counts'.
# alpha=2, beta=4 encodes a distribution with a weak peak at 1/3.
# alpha=4, beta=8 encodes a distribution with a stronger peak at 1/3.
# alpha=20, beta=40 encodes a distribution with a very strong peak at 1/3.

# Alpha and beta both need to be strictly positive, so we either have to sample them from a positive distribution, or apply a mapping.

# A model where alpha and beta are both linear functions of the covariates, and the intercepts for both alpha and beta are random
alpha = 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=(10.0, 3.0)),
        sigma=make_prior(dist_name="Normal", dist_params=(0.0, 3.0), mapping="softplus", mapping_params=(0.0, 3.0)),
    ),
    mapping="softplus",
    mapping_params=(0.0, 3.0),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
)
beta = 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=(10.0, 3.0)),
        sigma=make_prior(dist_name="Normal", dist_params=(0.0, 3.0), mapping="softplus", mapping_params=(0.0, 3.0)),
    ),
    mapping="softplus",
    mapping_params=(0.0, 3.0),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
)


beta_likelihood = BetaLikelihood(alpha, beta)

beta_HBR_model = HBR(
    name="beta_HBR_model",
    likelihood=beta_likelihood,
    cores=16,
    progressbar=True,
    draws=1500,
    tune=500,
    chains=4,
)

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=beta_HBR_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_beta/save_dir",
    # The scaler to use for the input data. Can be either one of "standardize", "minmax", "robminmax", "none"
    inscaler="minmax",
    # The scaler to use for the output data. Can be either one of "standardize", "minmax", "robminmax", "none"
    outscaler="minmax",
)

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 19 seconds

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.16	31
	2000	0	0.16	63
	2000	0	0.15	127
	2000	0	0.16	127

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 16 seconds

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.14	95
	2000	0	0.15	31
	2000	0	0.15	31
	2000	0	0.16	31

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 22 seconds

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.14	63
	2000	0	0.14	63
	2000	0	0.15	63
	2000	0	0.14	63

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 22 seconds

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	2000	0	0.14	63
	2000	0	0.15	63
	2000	0	0.15	63
	2000	0	0.14	127

<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.5741 ... -...
    centiles           (centile, observations, response_vars) float64 35kB 87...
    logp               (observations, response_vars) float64 7kB 1.161 ... 0.553
    Yhat               (observations, response_vars) float64 7kB 2.25e+03 ......
    statistics         (response_vars, statistic) float64 352B -0.4643 ... nan
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.5741 0.6577 ... -1.142 -1.132

    array([[ 5.74095715e-01,  6.57669282e-01, -8.74156542e-01,
             3.11343258e-01],
           [ 1.92828295e-01,  9.27751872e-01, -1.18763796e+00,
            -6.90106915e-01],
           [ 3.50756022e-01,  4.66806431e-01,  6.27048882e-01,
             3.06141763e-01],
           [ 3.13142547e-01,  1.28306414e+00, -8.27209895e-01,
             1.14139650e+00],
           [-1.37509618e+00, -9.06785218e-01, -8.06524893e-01,
                       -inf],
           [-7.32346384e-01, -4.13999399e-01,  1.26906938e-01,
            -3.82970154e-01],
           [ 5.80605763e-01,  2.65877825e+00,  2.96523907e-01,
             3.08793167e-01],
           [ 6.29676054e-01, -4.26208872e-01,  2.18605825e+00,
            -6.85277928e-01],
           [-9.92309055e-01,  7.48782299e-01, -9.90477127e-01,
            -6.82247800e-01],
           [-1.68418887e+00, -2.85582002e-01, -3.08493034e-01,
             3.57187231e-01],
    ...
           [-3.15872162e-01,  4.63500109e-01,  5.28420442e-01,
            -4.55544017e-01],
           [-1.48256116e+00, -4.51049474e-02, -4.19197747e-01,
            -5.54696582e-01],
           [-3.90789146e-01,  2.12656446e-01,  1.02566128e+00,
             6.04518578e-01],
           [ 2.40626556e-01,  9.17839756e-02, -1.75339393e-01,
            -7.89761839e-01],
           [-4.31455578e-01,  1.02369726e+00, -4.08275884e-01,
            -6.32281393e-01],
           [-2.36230831e+00, -4.00056018e-01,  3.99888468e-01,
             1.44207087e+00],
           [-1.55015157e+00,  3.21521635e-01, -9.59209289e-01,
            -8.66553270e-01],
           [-3.72074206e-01, -2.00771585e-01,  1.19584509e+00,
             3.25164696e-01],
           [ 3.17103069e+00,  1.11493763e-01,  2.26237369e+00,
             3.06359976e+00],
           [-1.06650535e+00,  1.12366856e+00, -1.14175552e+00,
            -1.13238397e+00]])

  • centiles
    (centile, observations, response_vars)
    float64
    879.8 2.911e+03 ... 6.107e+05

    array([[[8.79808160e+02, 2.91145804e+03, 1.33661101e+03, 3.98819762e+05],
            [5.84597691e+02, 3.14182955e+03, 1.55293556e+03, 4.75340556e+05],
            [4.29997133e+02, 3.09992048e+03, 1.46510690e+03, 4.33633970e+05],
            ...,
            [1.08499996e+03, 3.76785858e+03, 1.32078358e+03, 4.00200035e+05],
            [8.35420866e+02, 2.63710055e+03, 1.55778776e+03, 4.34024812e+05],
            [5.65283199e+02, 3.15324815e+03, 1.58018857e+03, 4.83504194e+05]],
    
           [[1.49815663e+03, 5.34178813e+03, 1.46640626e+03, 4.24689780e+05],
            [8.36328344e+02, 4.79232853e+03, 1.76141273e+03, 5.27977828e+05],
            [6.28807276e+02, 4.56205061e+03, 1.61879780e+03, 4.62683940e+05],
            ...,
            [1.98875772e+03, 1.11594498e+04, 1.46899506e+03, 4.37482012e+05],
            [1.09595347e+03, 4.00465456e+03, 1.76015529e+03, 4.67674630e+05],
            [8.19456865e+02, 4.80316607e+03, 1.76618525e+03, 5.22154842e+05]],
    
           [[2.10488534e+03, 8.66459928e+03, 1.56507118e+03, 4.46630834e+05],
            [1.06870579e+03, 6.59910861e+03, 1.90980756e+03, 5.65198964e+05],
            [8.30917832e+02, 6.13649086e+03, 1.73154411e+03, 4.85004652e+05],
            ...,
            [2.88393327e+03, 2.05101302e+04, 1.58332026e+03, 4.70625784e+05],
            [1.31364148e+03, 5.74231058e+03, 1.90409339e+03, 4.93666252e+05],
            [1.05843611e+03, 6.60730050e+03, 1.89791543e+03, 5.49496017e+05]],
    
           [[2.85081313e+03, 1.32368298e+04, 1.66832827e+03, 4.71167947e+05],
            [1.35039567e+03, 8.94612144e+03, 2.05336362e+03, 5.99832030e+05],
            [1.09078520e+03, 8.17748894e+03, 1.84576880e+03, 5.08398497e+05],
            ...,
            [3.96020185e+03, 2.84013143e+04, 1.70322058e+03, 5.07356506e+05],
            [1.56160691e+03, 8.18049990e+03, 2.04374089e+03, 5.20751406e+05],
            [1.35130244e+03, 8.94954192e+03, 2.02618226e+03, 5.75994548e+05]],
    
           [[4.11689516e+03, 2.06338009e+04, 1.82020175e+03, 5.09269746e+05],
            [1.84141557e+03, 1.30763636e+04, 2.23923026e+03, 6.40836368e+05],
            [1.56814518e+03, 1.18046892e+04, 2.00652576e+03, 5.42566222e+05],
            ...,
            [5.67076960e+03, 3.33442654e+04, 1.87795693e+03, 5.60607534e+05],
            [1.97024017e+03, 1.26996697e+04, 2.22596863e+03, 5.59564969e+05],
            [1.86633693e+03, 1.30682435e+04, 2.19581792e+03, 6.10686845e+05]]],
          shape=(5, 216, 4))

  • logp
    (observations, response_vars)
    float64
    1.161 0.3069 1.162 ... 0.5843 0.553

    array([[ 1.16109908e+00,  3.06927535e-01,  1.16184896e+00,
             1.22930628e+00],
           [ 2.41095075e+00,  7.18360940e-01,  4.42746454e-01,
             6.45821997e-01],
           [ 2.38976792e+00,  1.28726461e+00,  1.17624676e+00,
             1.29312353e+00],
           [ 1.80792705e+00, -3.60489911e-02, -2.25870880e-01,
             4.67836344e-01],
           [ 1.10762452e+00,  1.07087590e+00,  1.26881670e+00,
                       -inf],
           [ 2.44735801e+00,  1.63287959e+00,  1.06817573e+00,
             1.18957305e+00],
           [ 2.14203409e+00, -2.54794036e+00,  1.17804510e+00,
             1.15165506e+00],
           [ 2.07625223e+00,  1.58516846e+00, -1.02081328e+00,
             1.13373070e+00],
           [ 2.41287253e+00,  1.19551697e+00,  9.14191994e-01,
             1.18877352e+00],
           [ 1.55703751e+00,  1.43792281e+00,  1.12846796e+00,
             1.17605438e+00],
    ...
           [ 2.19915399e+00,  1.38000137e+00,  8.13448288e-01,
             4.63003680e-01],
           [ 1.71804332e+00,  1.54880619e+00,  1.27519839e+00,
             1.21398145e+00],
           [ 2.24754484e+00,  1.70844744e+00,  6.93333339e-01,
             6.44600821e-01],
           [ 2.62911052e+00,  1.47662005e+00,  1.34370982e+00,
             1.03085611e+00],
           [ 2.53092208e+00,  7.45078562e-01,  9.87134968e-01,
             1.01527378e+00],
           [ 7.25514287e-01,  1.79225542e+00,  1.20973614e+00,
             3.33343437e-01],
           [ 1.81406638e+00,  1.27600083e+00,  6.77683784e-01,
             4.94114876e-01],
           [ 1.15019949e+00, -2.06884226e-01,  5.58297003e-01,
             8.18467201e-01],
           [-2.89429236e+00,  1.36984744e+00, -1.05621462e+00,
            -3.26857735e+00],
           [ 2.21931033e+00,  4.64005696e-01,  5.84306452e-01,
             5.53017381e-01]])

  • Yhat
    (observations, response_vars)
    float64
    2.25e+03 9.814e+03 ... 5.486e+05

    array([[2.24967052e+03, 9.81372963e+03, 1.56998674e+03, 4.49359966e+05],
           [1.12194121e+03, 7.15437654e+03, 1.90475803e+03, 5.62586245e+05],
           [8.92928960e+02, 6.62056822e+03, 1.73311819e+03, 4.86144976e+05],
           [1.38171881e+03, 8.52937404e+03, 2.20280593e+03, 5.32870878e+05],
           [2.21172462e+03, 9.18129921e+03, 1.44483198e+03, 3.78583459e+05],
           [1.16686835e+03, 6.61898120e+03, 1.87417604e+03, 5.25639688e+05],
           [1.03316872e+03, 7.47954657e+03, 1.95855646e+03, 5.38120875e+05],
           [5.26078758e+02, 6.71307968e+03, 1.62546272e+03, 4.92364763e+05],
           [1.03021359e+03, 5.77663510e+03, 1.63530390e+03, 4.84538822e+05],
           [1.27909138e+03, 7.27014587e+03, 1.54920211e+03, 4.34977534e+05],
           [1.16929411e+03, 6.39934212e+03, 1.87339950e+03, 5.39823802e+05],
           [9.86187648e+02, 6.30822159e+03, 1.57238672e+03, 4.62717705e+05],
           [1.03176109e+03, 5.72668907e+03, 1.63364741e+03, 4.86532626e+05],
           [1.24363516e+03, 5.67561086e+03, 1.67254168e+03, 4.50090250e+05],
           [1.21690783e+03, 8.42661530e+03, 2.09722435e+03, 5.30958245e+05],
           [9.01593095e+02, 6.33973486e+03, 1.66823447e+03, 5.06633733e+05],
           [1.45801540e+03, 7.47274308e+03, 1.83582716e+03, 5.48062800e+05],
           [1.03422958e+03, 5.66601517e+03, 1.63011292e+03, 4.88515648e+05],
           [1.20892663e+03, 7.71907283e+03, 2.01603094e+03, 4.91613619e+05],
           [1.24327655e+03, 5.21563473e+03, 1.66927923e+03, 4.57508912e+05],
    ...
           [7.60321942e+02, 6.29665047e+03, 1.67771273e+03, 5.09613863e+05],
           [1.34891229e+03, 6.13879515e+03, 1.89899778e+03, 5.05432095e+05],
           [1.29416778e+03, 6.38577418e+03, 1.61375193e+03, 4.94531095e+05],
           [1.71118885e+03, 6.48559736e+03, 1.54541241e+03, 4.27783424e+05],
           [1.19207424e+03, 7.99539745e+03, 1.87883644e+03, 5.21951269e+05],
           [8.43516079e+03, 2.82477424e+04, 1.61808669e+03, 4.88361002e+05],
           [1.23577515e+03, 5.49241934e+03, 1.67819392e+03, 4.54636444e+05],
           [6.37047239e+02, 8.01482547e+03, 1.85228473e+03, 5.39113188e+05],
           [1.15463605e+03, 7.51132814e+03, 1.93854528e+03, 5.33816042e+05],
           [1.03422958e+03, 5.66601517e+03, 1.63011292e+03, 4.88515648e+05],
           [1.55092109e+03, 6.21104942e+03, 1.62765935e+03, 5.16710431e+05],
           [1.29416778e+03, 6.38577418e+03, 1.61375193e+03, 4.94531095e+05],
           [1.38146571e+03, 5.54452374e+03, 1.44377487e+03, 4.30296805e+05],
           [5.20957728e+02, 6.61813182e+03, 1.63064721e+03, 4.97159505e+05],
           [1.16727786e+03, 6.48046211e+03, 1.87398145e+03, 5.34612095e+05],
           [1.03191177e+03, 5.88288098e+03, 1.63242761e+03, 4.78871278e+05],
           [1.12140454e+03, 7.11368556e+03, 1.90447928e+03, 5.66428715e+05],
           [3.06564266e+03, 1.96433860e+04, 1.58923207e+03, 4.74220845e+05],
           [1.34656848e+03, 6.45054606e+03, 1.89959817e+03, 4.94818015e+05],
           [1.11648384e+03, 7.16017187e+03, 1.89426617e+03, 5.48612541e+05]])

  • statistics
    (response_vars, statistic)
    float64
    -0.4643 0.02093 ... 0.6457 nan

    array([[-4.64267626e-01,  2.09259259e-02,  3.84606326e-01,
            -9.74018094e+00, -1.91563079e+00, -4.95971861e-01,
             7.40258648e+02,  4.67599655e-01,  3.92603185e-13,
             1.49597186e+00,  9.93927473e-01],
           [ 1.10962109e-01,  1.94444444e-02,  4.31789035e-01,
            -1.06826277e+01, -9.92121932e-01,  1.10107299e-01,
             3.68947705e+03,  2.69322160e-01,  6.08264063e-05,
             8.89892701e-01,  9.78975661e-01],
           [ 2.45820617e-01,  2.50000000e-02,  9.73796170e-02,
            -7.53263171e+00, -6.50105676e-01,  2.42966527e-01,
             2.05287944e+02,  4.99758847e-01,  4.76553315e-15,
             7.57033473e-01,  9.90262088e-01],
           [ 3.63155733e-01,  2.59259259e-02,  6.50527730e-02,
                        inf,             inf,  3.54288779e-01,
             3.93369730e+04,  6.34723158e-01,  9.32846749e-26,
             6.45711221e-01,             nan]])

• 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 0x141ffdf80>, {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 2. plot_qq: Plot the QQ-plot of the predicted Z-scores 3. 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.
)

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.478848	0.013039	0.059172	-13.107157	-0.826584	0.478654	37645.305005	0.703202	1.547316e-129	0.521346	0.983431
Right-Amygdala	0.385791	0.022738	0.088500	-7.693682	-0.773487	0.385457	192.062106	0.604786	4.163637e-87	0.614543	0.977994
Right-Lateral-Ventricle	0.165282	0.026032	0.395923	-10.756175	-1.081352	0.161034	3526.609410	0.393598	2.528902e-33	0.838966	0.921989
WM-hypointensities	0.011507	0.030905	0.343233	-10.110444	-1.986224	0.003670	815.207354	0.515763	9.270689e-60	0.996330	0.971463

statistic	EXPV	MACE	MAPE	MSLL	NLL	R2	RMSE	Rho	Rho_p	SMSE	ShapiroW
response_vars
CortexVol	0.363156	0.025926	0.065053	inf	inf	0.354289	39336.972955	0.634723	9.328467e-26	0.645711	NaN
Right-Amygdala	0.245821	0.025000	0.097380	-7.532632	-0.650106	0.242967	205.287944	0.499759	4.765533e-15	0.757033	0.990262
Right-Lateral-Ventricle	0.110962	0.019444	0.431789	-10.682628	-0.992122	0.110107	3689.477050	0.269322	6.082641e-05	0.889893	0.978976
WM-hypointensities	-0.464268	0.020926	0.384606	-9.740181	-1.915631	-0.495972	740.258648	0.467600	3.926032e-13	1.495972	0.993927

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.5741 ... -...
    logp               (observations, response_vars) float64 7kB 1.161 ... 0.553
    Yhat               (observations, response_vars) float64 7kB 2.25e+03 ......
    statistics         (response_vars, statistic) float64 352B -0.4643 ... nan
    centiles           (centile, observations, response_vars) float64 35kB 87...
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.5741 0.6577 ... -1.142 -1.132

    array([[ 5.74095715e-01,  6.57669282e-01, -8.74156542e-01,
             3.11343258e-01],
           [ 1.92828295e-01,  9.27751872e-01, -1.18763796e+00,
            -6.90106915e-01],
           [ 3.50756022e-01,  4.66806431e-01,  6.27048882e-01,
             3.06141763e-01],
           [ 3.13142547e-01,  1.28306414e+00, -8.27209895e-01,
             1.14139650e+00],
           [-1.37509618e+00, -9.06785218e-01, -8.06524893e-01,
                       -inf],
           [-7.32346384e-01, -4.13999399e-01,  1.26906938e-01,
            -3.82970154e-01],
           [ 5.80605763e-01,  2.65877825e+00,  2.96523907e-01,
             3.08793167e-01],
           [ 6.29676054e-01, -4.26208872e-01,  2.18605825e+00,
            -6.85277928e-01],
           [-9.92309055e-01,  7.48782299e-01, -9.90477127e-01,
            -6.82247800e-01],
           [-1.68418887e+00, -2.85582002e-01, -3.08493034e-01,
             3.57187231e-01],
    ...
           [-3.15872162e-01,  4.63500109e-01,  5.28420442e-01,
            -4.55544017e-01],
           [-1.48256116e+00, -4.51049474e-02, -4.19197747e-01,
            -5.54696582e-01],
           [-3.90789146e-01,  2.12656446e-01,  1.02566128e+00,
             6.04518578e-01],
           [ 2.40626556e-01,  9.17839756e-02, -1.75339393e-01,
            -7.89761839e-01],
           [-4.31455578e-01,  1.02369726e+00, -4.08275884e-01,
            -6.32281393e-01],
           [-2.36230831e+00, -4.00056018e-01,  3.99888468e-01,
             1.44207087e+00],
           [-1.55015157e+00,  3.21521635e-01, -9.59209289e-01,
            -8.66553270e-01],
           [-3.72074206e-01, -2.00771585e-01,  1.19584509e+00,
             3.25164696e-01],
           [ 3.17103069e+00,  1.11493763e-01,  2.26237369e+00,
             3.06359976e+00],
           [-1.06650535e+00,  1.12366856e+00, -1.14175552e+00,
            -1.13238397e+00]])

  • logp
    (observations, response_vars)
    float64
    1.161 0.3069 1.162 ... 0.5843 0.553

    array([[ 1.16109908e+00,  3.06927535e-01,  1.16184896e+00,
             1.22930628e+00],
           [ 2.41095075e+00,  7.18360940e-01,  4.42746454e-01,
             6.45821997e-01],
           [ 2.38976792e+00,  1.28726461e+00,  1.17624676e+00,
             1.29312353e+00],
           [ 1.80792705e+00, -3.60489911e-02, -2.25870880e-01,
             4.67836344e-01],
           [ 1.10762452e+00,  1.07087590e+00,  1.26881670e+00,
                       -inf],
           [ 2.44735801e+00,  1.63287959e+00,  1.06817573e+00,
             1.18957305e+00],
           [ 2.14203409e+00, -2.54794036e+00,  1.17804510e+00,
             1.15165506e+00],
           [ 2.07625223e+00,  1.58516846e+00, -1.02081328e+00,
             1.13373070e+00],
           [ 2.41287253e+00,  1.19551697e+00,  9.14191994e-01,
             1.18877352e+00],
           [ 1.55703751e+00,  1.43792281e+00,  1.12846796e+00,
             1.17605438e+00],
    ...
           [ 2.19915399e+00,  1.38000137e+00,  8.13448288e-01,
             4.63003680e-01],
           [ 1.71804332e+00,  1.54880619e+00,  1.27519839e+00,
             1.21398145e+00],
           [ 2.24754484e+00,  1.70844744e+00,  6.93333339e-01,
             6.44600821e-01],
           [ 2.62911052e+00,  1.47662005e+00,  1.34370982e+00,
             1.03085611e+00],
           [ 2.53092208e+00,  7.45078562e-01,  9.87134968e-01,
             1.01527378e+00],
           [ 7.25514287e-01,  1.79225542e+00,  1.20973614e+00,
             3.33343437e-01],
           [ 1.81406638e+00,  1.27600083e+00,  6.77683784e-01,
             4.94114876e-01],
           [ 1.15019949e+00, -2.06884226e-01,  5.58297003e-01,
             8.18467201e-01],
           [-2.89429236e+00,  1.36984744e+00, -1.05621462e+00,
            -3.26857735e+00],
           [ 2.21931033e+00,  4.64005696e-01,  5.84306452e-01,
             5.53017381e-01]])

  • Yhat
    (observations, response_vars)
    float64
    2.25e+03 9.814e+03 ... 5.486e+05

    array([[2.24967052e+03, 9.81372963e+03, 1.56998674e+03, 4.49359966e+05],
           [1.12194121e+03, 7.15437654e+03, 1.90475803e+03, 5.62586245e+05],
           [8.92928960e+02, 6.62056822e+03, 1.73311819e+03, 4.86144976e+05],
           [1.38171881e+03, 8.52937404e+03, 2.20280593e+03, 5.32870878e+05],
           [2.21172462e+03, 9.18129921e+03, 1.44483198e+03, 3.78583459e+05],
           [1.16686835e+03, 6.61898120e+03, 1.87417604e+03, 5.25639688e+05],
           [1.03316872e+03, 7.47954657e+03, 1.95855646e+03, 5.38120875e+05],
           [5.26078758e+02, 6.71307968e+03, 1.62546272e+03, 4.92364763e+05],
           [1.03021359e+03, 5.77663510e+03, 1.63530390e+03, 4.84538822e+05],
           [1.27909138e+03, 7.27014587e+03, 1.54920211e+03, 4.34977534e+05],
           [1.16929411e+03, 6.39934212e+03, 1.87339950e+03, 5.39823802e+05],
           [9.86187648e+02, 6.30822159e+03, 1.57238672e+03, 4.62717705e+05],
           [1.03176109e+03, 5.72668907e+03, 1.63364741e+03, 4.86532626e+05],
           [1.24363516e+03, 5.67561086e+03, 1.67254168e+03, 4.50090250e+05],
           [1.21690783e+03, 8.42661530e+03, 2.09722435e+03, 5.30958245e+05],
           [9.01593095e+02, 6.33973486e+03, 1.66823447e+03, 5.06633733e+05],
           [1.45801540e+03, 7.47274308e+03, 1.83582716e+03, 5.48062800e+05],
           [1.03422958e+03, 5.66601517e+03, 1.63011292e+03, 4.88515648e+05],
           [1.20892663e+03, 7.71907283e+03, 2.01603094e+03, 4.91613619e+05],
           [1.24327655e+03, 5.21563473e+03, 1.66927923e+03, 4.57508912e+05],
    ...
           [7.60321942e+02, 6.29665047e+03, 1.67771273e+03, 5.09613863e+05],
           [1.34891229e+03, 6.13879515e+03, 1.89899778e+03, 5.05432095e+05],
           [1.29416778e+03, 6.38577418e+03, 1.61375193e+03, 4.94531095e+05],
           [1.71118885e+03, 6.48559736e+03, 1.54541241e+03, 4.27783424e+05],
           [1.19207424e+03, 7.99539745e+03, 1.87883644e+03, 5.21951269e+05],
           [8.43516079e+03, 2.82477424e+04, 1.61808669e+03, 4.88361002e+05],
           [1.23577515e+03, 5.49241934e+03, 1.67819392e+03, 4.54636444e+05],
           [6.37047239e+02, 8.01482547e+03, 1.85228473e+03, 5.39113188e+05],
           [1.15463605e+03, 7.51132814e+03, 1.93854528e+03, 5.33816042e+05],
           [1.03422958e+03, 5.66601517e+03, 1.63011292e+03, 4.88515648e+05],
           [1.55092109e+03, 6.21104942e+03, 1.62765935e+03, 5.16710431e+05],
           [1.29416778e+03, 6.38577418e+03, 1.61375193e+03, 4.94531095e+05],
           [1.38146571e+03, 5.54452374e+03, 1.44377487e+03, 4.30296805e+05],
           [5.20957728e+02, 6.61813182e+03, 1.63064721e+03, 4.97159505e+05],
           [1.16727786e+03, 6.48046211e+03, 1.87398145e+03, 5.34612095e+05],
           [1.03191177e+03, 5.88288098e+03, 1.63242761e+03, 4.78871278e+05],
           [1.12140454e+03, 7.11368556e+03, 1.90447928e+03, 5.66428715e+05],
           [3.06564266e+03, 1.96433860e+04, 1.58923207e+03, 4.74220845e+05],
           [1.34656848e+03, 6.45054606e+03, 1.89959817e+03, 4.94818015e+05],
           [1.11648384e+03, 7.16017187e+03, 1.89426617e+03, 5.48612541e+05]])

  • statistics
    (response_vars, statistic)
    float64
    -0.4643 0.02093 ... 0.6457 nan

    array([[-4.64267626e-01,  2.09259259e-02,  3.84606326e-01,
            -9.74018094e+00, -1.91563079e+00, -4.95971861e-01,
             7.40258648e+02,  4.67599655e-01,  3.92603185e-13,
             1.49597186e+00,  9.93927473e-01],
           [ 1.10962109e-01,  1.94444444e-02,  4.31789035e-01,
            -1.06826277e+01, -9.92121932e-01,  1.10107299e-01,
             3.68947705e+03,  2.69322160e-01,  6.08264063e-05,
             8.89892701e-01,  9.78975661e-01],
           [ 2.45820617e-01,  2.50000000e-02,  9.73796170e-02,
            -7.53263171e+00, -6.50105676e-01,  2.42966527e-01,
             2.05287944e+02,  4.99758847e-01,  4.76553315e-15,
             7.57033473e-01,  9.90262088e-01],
           [ 3.63155733e-01,  2.59259259e-02,  6.50527730e-02,
                        inf,             inf,  3.54288779e-01,
             3.93369730e+04,  6.34723158e-01,  9.32846749e-26,
             6.45711221e-01,             nan]])

  • centiles
    (centile, observations, response_vars)
    float64
    879.8 2.911e+03 ... 6.107e+05

    array([[[8.79808160e+02, 2.91145804e+03, 1.33661101e+03, 3.98819762e+05],
            [5.84597691e+02, 3.14182955e+03, 1.55293556e+03, 4.75340556e+05],
            [4.29997133e+02, 3.09992048e+03, 1.46510690e+03, 4.33633970e+05],
            ...,
            [1.08499996e+03, 3.76785858e+03, 1.32078358e+03, 4.00200035e+05],
            [8.35420866e+02, 2.63710055e+03, 1.55778776e+03, 4.34024812e+05],
            [5.65283199e+02, 3.15324815e+03, 1.58018857e+03, 4.83504194e+05]],
    
           [[1.49815663e+03, 5.34178813e+03, 1.46640626e+03, 4.24689780e+05],
            [8.36328344e+02, 4.79232853e+03, 1.76141273e+03, 5.27977828e+05],
            [6.28807276e+02, 4.56205061e+03, 1.61879780e+03, 4.62683940e+05],
            ...,
            [1.98875772e+03, 1.11594498e+04, 1.46899506e+03, 4.37482012e+05],
            [1.09595347e+03, 4.00465456e+03, 1.76015529e+03, 4.67674630e+05],
            [8.19456865e+02, 4.80316607e+03, 1.76618525e+03, 5.22154842e+05]],
    
           [[2.10488534e+03, 8.66459928e+03, 1.56507118e+03, 4.46630834e+05],
            [1.06870579e+03, 6.59910861e+03, 1.90980756e+03, 5.65198964e+05],
            [8.30917832e+02, 6.13649086e+03, 1.73154411e+03, 4.85004652e+05],
            ...,
            [2.88393327e+03, 2.05101302e+04, 1.58332026e+03, 4.70625784e+05],
            [1.31364148e+03, 5.74231058e+03, 1.90409339e+03, 4.93666252e+05],
            [1.05843611e+03, 6.60730050e+03, 1.89791543e+03, 5.49496017e+05]],
    
           [[2.85081313e+03, 1.32368298e+04, 1.66832827e+03, 4.71167947e+05],
            [1.35039567e+03, 8.94612144e+03, 2.05336362e+03, 5.99832030e+05],
            [1.09078520e+03, 8.17748894e+03, 1.84576880e+03, 5.08398497e+05],
            ...,
            [3.96020185e+03, 2.84013143e+04, 1.70322058e+03, 5.07356506e+05],
            [1.56160691e+03, 8.18049990e+03, 2.04374089e+03, 5.20751406e+05],
            [1.35130244e+03, 8.94954192e+03, 2.02618226e+03, 5.75994548e+05]],
    
           [[4.11689516e+03, 2.06338009e+04, 1.82020175e+03, 5.09269746e+05],
            [1.84141557e+03, 1.30763636e+04, 2.23923026e+03, 6.40836368e+05],
            [1.56814518e+03, 1.18046892e+04, 2.00652576e+03, 5.42566222e+05],
            ...,
            [5.67076960e+03, 3.33442654e+04, 1.87795693e+03, 5.60607534e+05],
            [1.97024017e+03, 1.26996697e+04, 2.22596863e+03, 5.59564969e+05],
            [1.86633693e+03, 1.30682435e+04, 2.19581792e+03, 6.10686845e+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 0x141ffdf80>, {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(test, reference_batch_effect=reference_batch_effect)  # <- easy

plt.style.use("seaborn-v0_8")
df = test.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,
    show_other_data=True,
    harmonize_data=True,
    show_legend=True,
)

# 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=True,  # Show data points that do not match any batch effects
    show_centile_labels=False,
    harmonize_data=True,  # Set this to False to see the difference
    show_legend=False,  # Don't show the legend because it crowds the plot
)

/opt/anaconda3/envs/uv_refactor/lib/python3.12/site-packages/pcntoolkit/util/output.py:218: UserWarning: Process: 64750 - 2025-06-23 16:05:58 - remove_Nan is set to False. Ensure your data does not contain NaNs in critical columns, or handle them appropriately.
  warnings.warn(message)

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 are some more HBR model configuration examples. You can replace the template_hbr that we used in the example with any of these.

# The beta likelihood is a bit different from the other likelihoods, it assumes the response variable is bounded between 0 and 1. The means that we have to apply a min max scaling to the response variable before fitting the model. Luckily, this is built into the PCNtoolkit, we just have to set 'outscaler' to 'minmax' in the `NormativeModel` constructor.

# The beta likelihood takes two parameters, alpha and beta. Intuitively, these encode 'pseudo-counts'.
# alpha=2, beta=4 encodes a distribution with a weak peak at 1/3.
# alpha=4, beta=8 encodes a distribution with a stronger peak at 1/3.
# alpha=20, beta=40 encodes a distribution with a very strong peak at 1/3.

# Alpha and beta both need to be strictly positive, so we either have to sample them from a positive distribution, or apply a mapping.


# Pick any of these three models for alpha and beta

# A very simple model where alpha and beta are both constant
alpha = make_prior(
    dist_name="Gamma",
    dist_params=(3.0, 1.0),
)
beta = make_prior(
    dist_name="Gamma",
    dist_params=(3.0, 1.0),
)

# A model where alpha and beta are both linear functions of the covariates
alpha = make_prior(
    linear=True,
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 10.0)),
    intercept=make_prior(dist_name="Normal", dist_params=(10.0, 3.0)),
    mapping="softplus",
    mapping_params=(0.0, 3.0),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
)
beta = make_prior(
    linear=True,
    slope=make_prior(dist_name="Normal", dist_params=(0.0, 10.0)),
    intercept=make_prior(dist_name="Normal", dist_params=(10.0, 3.0)),
    mapping="softplus",
    mapping_params=(0.0, 3.0),
    basis_function=BsplineBasisFunction(basis_column=0, nknots=5, degree=3),
)


beta_likelihood = BetaLikelihood(alpha, beta)

beta_HBR_model = HBR(
    name="beta_HBR_model",
    likelihood=beta_likelihood,
    cores=16,
    progressbar=True,
    draws=1500,
    tune=500,
    chains=4,
)