This paper introduces an extension of hierarchical Bayesian regression (HBR) to model non-Gaussian data with heteroskedastic skewness and kurtosis, using the sinh-arcsinh (SHASH) distribution. The authors address the limitations of traditional Gaussian assumptions in normative modelling, which can be restrictive in certain applications. They propose a novel reparameterization and Markov chain Monte Carlo (MCMC) sampling approach to perform inference in this model. Using a large neuroimaging dataset collected at 82 different sites, they demonstrate that the extended HBR framework with a SHASH likelihood outperforms or matches the performance of a warped Bayesian linear regression baseline model on most datasets. The method is shown to provide better control over the parameters governing the shape of distributions and is useful for accurately modeling highly nonlinear relationships between aging and imaging-derived phenotypes. The contributions of this work include a reparameterization of the SHASH distribution, a thorough analysis of the HBR method with a SHASH likelihood, extensive evaluation on a multi-site neuroimaging dataset, and an extension to the existing HBR implementation in the pcnt toolkit. The paper also discusses the advantages and limitations of the proposed approach and provides a detailed theoretical background, experimental results, and ethical considerations.This paper introduces an extension of hierarchical Bayesian regression (HBR) to model non-Gaussian data with heteroskedastic skewness and kurtosis, using the sinh-arcsinh (SHASH) distribution. The authors address the limitations of traditional Gaussian assumptions in normative modelling, which can be restrictive in certain applications. They propose a novel reparameterization and Markov chain Monte Carlo (MCMC) sampling approach to perform inference in this model. Using a large neuroimaging dataset collected at 82 different sites, they demonstrate that the extended HBR framework with a SHASH likelihood outperforms or matches the performance of a warped Bayesian linear regression baseline model on most datasets. The method is shown to provide better control over the parameters governing the shape of distributions and is useful for accurately modeling highly nonlinear relationships between aging and imaging-derived phenotypes. The contributions of this work include a reparameterization of the SHASH distribution, a thorough analysis of the HBR method with a SHASH likelihood, extensive evaluation on a multi-site neuroimaging dataset, and an extension to the existing HBR implementation in the pcnt toolkit. The paper also discusses the advantages and limitations of the proposed approach and provides a detailed theoretical background, experimental results, and ethical considerations.