x. D: GAMLSS allows to model y as a multiparameter non-gaussian distribution, where location, scale, and shape parameters can be modeled as functions of input variables.

A linear model can be described as = β + ϵ, ϵ ∼ (0, σ 2 ) or using a slightly different notation ∼ (µ, σ 2 ), µ = β where y is the target variable, βis a vector of estimated coefficients, X is a design matrix and ϵ is an error vector, i.e., y depends linearly on the predictors X and the error term is assumed to be normally distributed. GAM extends LM by relaxing the linearity assumption: ∼ (µ, σ 2 ), µ = 1( 1) + 2( 2) +... + ( )

Where ( ) is a possibly non-linear or non-parametric function of . The original formulation of GAMs (Hastie and Tibshirani 1990) allowed for almost arbitrary functions, including loess smoothing, decision trees, or neural networks that could then be an additive part of a larger model. More modern interpretations of GAMs (Ref: Wood) limit these functions to smooth mostly spline-based functions of X, where the complexity of the model can be controlled by penalizing the wiggliness of a function (In GAM nomenclature, ‘wiggliness’ is given a precise technical meaning usually referring to an integrated second derivative of a function). The wiggliness of the function is a parameter that can be estimated from the data, which allows for flexible modeling of non-linear effects with an automatic selection of model complexity. In a linear model, we can model non-linear effects by a spline basis expansion of the input variables, and the number of knots will control the model complexity. In GAM, the number of knots is usually fixed at some high number, and the model complexity is controlled by directly penalizing the second derivative (i.e., the wiggliness) of the function, which can be implemented as a quadratic penalty on weights of the spline basis expansion.

GAMLSS is an extension of GAM introduced by (Stasinopoulos et al. 2017; Rigby and Stasinopoulos 2005) , which allows to model multiple parameters of a multiple parameter distributions as additive functions, i.e., instead of just modeling the expected mean of the target variable and assuming that the distribution around this mean is normally distributed, we can model location, scale, and shape parameters of the distribution. So for a four-parameter distribution ∼ (µ, σ, ν, τ) µ = µ,1( 1) + µ,2( 2)+ …+ µ, ( )

σ = σ,1( 1) + σ,2( 2)+ …+ σ, ( )

ν = ν,1( 1) + ν,2( 2)+ …+ ν, ( )

τ = τ,1( 1) + τ,2( 2)+ …+ τ, ( )

Here D is a parametric distribution with parameters µ, σ, ν, τ, and each of the parameters can be modeled as a separate function of input variables. This allows for accounting for heteroskedasticity by directly modeling the scale of a distribution (σ) as a function of other variables or modeling shapes of non-gaussian distributions. Each parameter can be modeled based on a different set of variables and using a completely different structural form. So it is possible to have a complex model for µ but at the same time a very restricted model for σ.

GAMLSS is usually estimated using a penalized maximum likelihood estimation. The smoothing functions () can be expressed as γ where is a matrix of basis expansion of , and γ is a set of coefficients penalized using a quadratic penaltyλγ γ where λ is a smoothing hyperparameter. We can also say thatγare random effects, assumed to be distributed as γ ∼ (0, −). If we are estimating a model for 4 parameter distribution, that depends on J smooth functions, the criterion for the penalized maximum likelihood is = − where = ∑ =1

4 12 ∑ =1 =∑1 λ γ

γ ( | µ , σ , ν , τ )

Traditionally, this is done iteratively using a Newton-Raphson or Fisher scoring (Rigby and Stasinopoulos 2005) , with recent algorithmic improvements improving computational complexity of these models and allowing their use in big datasets (Wood and Fasiolo 2017) . Full Bayesian inference using a posterior sampling is also possible (Umlauf, Klein, and Zeileis 2018) .

In the rest of this paper, we are using SHASH (Jones and Pewsey 2009) distribution as implemented in the mgcv library (Wood 2017) as an example multi-parameter distribution that can be fitted using the GAMLSS approach. The four parameters roughly correspond to location, scale, skewness, and kurtosis of the distribution. SHASH has the normal distribution as its special cases, but based on its parameters, it can also be skewed or more or less kurtic than the normal distribution. where , ,and Parameters µ, σ, ν, τcontrol location, scale, skewness, and kurtosis respectively.

The most important measures are the measures of total goodness of fit that are sensitive to all possible misspecification of the model, including badly estimated scale or shape of the distribution, heteroskedasticity, or differences between linear or nonlinear models. For this reason, various popular measures used to validate regression models, such as mean absolute error (MAE), R2, correlation, or root mean squared error (RMSE), are not appropriate since these are calculated based on the predictive meanbut do not take into account other aspects of modeled distribution. Appropriate measures are so-called proper scoring rules (Gneiting, Balabdaoui, and Raftery 2005) , which are maximized by the true data-generating model. One possible measure is the logarithmic score or validation global deviance (Stasinopoulos et al. 2017) or sometimes predictive deviance. =− 2() =− 2 (|) where is the model likelihood, and (|) is a predicted probability density by the model, at the target variable y in the test set. It is called a “validation global deviance”, since the model is estimated in a training set, and the deviance is calculated based on the data in the test set. The logarithmic score is defined similarly = () = (|) which can be scaled to aid interpretation as for example, deviance explained or mean standardized log loss (MSLL) (Rasmussen and Williams 2005) , which is averaged across all data points. The MSLL is computed as follows =− () − () =− (|) − (|) Where is a reference model, usually a trivial model that predicts only an mean and standard deviation of the sample. We can see that these measures are just linear transformations of each other, therefore the relative ranking of models will be the same, provided the reference model used to calculate MSLL will not change. The deviance is the predicted density of an individual's observed value, and therefore these measures depend on the location, scale, and shape of the predicted distribution.

One potential drawback is that these measures are sensitive to outliers, and unless care is taken with regard to treatment of outliers, this might lead to selecting the wrong model. For example, removal or Winsorizing of outliers will lead to models with underestimated variance and kurtosis of the distribution. One way to make these measures more resistant to outliers while keeping them proper is to use a censored likelihood (Diks, Panchenko, and van Dijk 2011) , which modifies likelihood definition to where where C is a chosen censored region. So if a predicted density is not in the censored region, the likelihood is calculated, as usual, otherwise, the log probability of a y being anywhere in the censored region given the model is used instead. This censored likelihood can be used instead of standard likelihood in other likelihood based measures (e.g., VGD or MSLL mentioned above), to create their robust versions.

For example, we might decide to censor values that are not between the 1st and 99th percentile. If the target value y for a subject is anywhere between the 1st and 99th percentile, the usual logarithmic score will be calculated. If the y is anywhere else, regardless of if the subject is an extreme outlier or not, the resulting score will be log(0.02), since the probability of y being outside of the 1st-99th percentile range is 2 percent. This means that this score would prefer models where a shape of the distribution between the 1st-99th percentiles is modeled accurately, and 2% of subjects are outside of this interval, but the shape of the distribution outside of the interval would not matter. The censored likelihood can also be used if we want to evaluate if some specific part of the distribution is modeled correctly, for example, only the distribution's right tail. For example, if we want to detect a disorder, and we know that the disorder has a neurodegenerative effect, we might want to focus on accurate modeling of the right tail of the distribution since this is where we expect clinical cases to lie. Although likelihood based measures are affected by all aspects of model fitting, it is still worth validating specific aspects of model fitting such as calibration and central tendency separately.

Measures of goodness of fit of the central tendency measure how well the center of the distribution is predicted. This is mainly a concern when comparing which model is better at modeling a possibly non-linear trend in the data. If the central tendency is not fitted correctly, we might expect, for example, that z-scores for younger and older subjects are systematically underestimated, and middle-aged subjects are overestimated. The center of the distribution can refer to the mean, the median, or the mode of the predicted distribution, and different interpretations may be preferred, for example, depending on whether the distribution is symmetric. The measures include the traditional measures commonly used to evaluate the performance of regression models, such as MAE, correlation, R2, or RMSE. There is nothing inherently better or worse about these measures; they all have their place, however, there are a few potential pitfalls when evaluating non-Gaussian models. In a Gaussian model, mean, median, and mode are the same, and they are fitted to minimize MSE. This means that MSE will often prefer Gaussian models, even if non-gaussian models perform better. The safest option is to use the predicted median and evaluate it using MAE, which is minimized by the median (Hanley et al. 2001) . The predicted mean might be technically hard to compute since if the mean is not one of the estimated parameters of a distribution, we would have to numerically integrate the whole density function to calculate it. Measures based on mean might also be negatively affected by skewness or outliers in the data.

Calibration assesses if the estimated centiles of a distribution match the observed centiles, or if the modeled distribution's shape is well approximating the observed distribution. Calibration in normative models is best evaluated using quantile randomized residuals, also known as randomized z-scores or pseudo-z-scores (Umlauf, Klein, and Zeileis 2018; Dunn and Smyth 1996) . These are quantiles of the fitted distribution mapped to z-scores of a standard normal distribution corresponding to the same quantiles. So if a subject is in a 95th percentile according to the GAMLSS model, their randomized z-score will be 1.97 (because 95 percentile of a standard normal distribution is at 1.97), regardless of the shape of the modeled distribution.

If the modeled distribution matches the observed distribution well, the randomized z-scores should be by construction normally distributed, regardless of the modeled distribution's shape. The most informative are visual displays of the residuals in the form of either histogram, Q-Q plots, P-P plots, or wormplots (detrended Q-Q plot) (Hanley et al. 2001; van Buuren and Fredriks 2001) . Various ways to visualize deviations from the assumed distribution are shown in figure 2. Of course, visual inspection is not always possible, for example, when modeling hundreds or even hundreds of thousands of voxels, therefore we need to have quantitative summary measures as well. Since the pseudo-z-scores should be normally distributed, calibration can be measured as various forms of deviations from the normality of the distribution pseudo-z-scores. One possibility is to use the test statistics from the Shapiro-Wilk test of normality (Shapiro and Wilk 1965) or a test statistic from other normality tests. In this case, we are not interested in calculating a p-value from the test but in using the test statistics to quantify the degree to which the observed distribution deviates from the assumed normal distribution. W is usually close to 1. Perfectly normally distributed data would obtain the value of W of 1, while the values less than 1 indicate deviation from normality. We can also measure scale, skewness, and kurtosis of the pseudo-z-scores separately to give us more insight into how the model is potentially misspecified.

Heteroskedasticity means that a variance of the response variable changes as a function of some variable in the dataset, for example, with age or scan-site. If this aspect is not modeled correctly, some subjects' abnormality will be systematically overestimated, and others underestimated. To investigate heteroskedasticity, we can again examine randomized residuals, which should by construction have a standard deviation of 1, the standard deviation should not be related to any variables in the data, i.e., in each scan-site or age bin. There are many ways how heteroskedasticity can be detected. One possibility is to examine the data visually by plotting squared residuals versus predicted variables and fitting a curve to this plot. If the residuals are homoskedastic, this curve should be close to a straight line. We can also calculate the goodness of fit of this curve to use as a quantitative measure of heteroskedasticity.

In the following section, we will use data from UK Biobank UKB (Miller et al. 2016) . We used regional cortical thickness data obtained using a processing pipeline as described elsewhere (Alfaro-Almagro et al. 2018) . For simplicity, we will start by modeling the effect of age on median cortical thickness in women. We split data into 70% training set and 30% test set. If this was a smaller dataset, we might consider using cross-validation instead of one test set. We start by fitting several models: a Gaussian model with linear effects (M1a), a non-linear Gaussian homoskedastic model (M1b), a non-linear Gaussian heteroskedastic model (M1c), and a non-linear heteroskedastic SHASH model made available under aCC-BY 4.0 International license. (non-Gaussian)(M1d), and a non-linear heteroskedastic SHASH model where location, scale, and shape also depends on age (Mm1e),: M1a: ℎ ∼ (µ, σ), = β M1b: ℎ ∼ (µ, σ), = M1c: ℎ ∼ (µ, σ), = M1d: ℎ ∼ M1e: ℎ ∼ () It is recommended to consider models with increased complexity in the order specified above, and not, for example, to consider models that have more complicated models for shape parameters than they have for scale or location (Stasinopoulos et al. 2017) . Therefore, we are not going to fit a model with for example a linear effect of age on µ but non-linear on σ. Figure 4 shows resulting model fits, and table 1 contains test set deviance and training set AIC and BIC. We can see that the effect of age seems to be slightly non-linear, while the distribution is slightly heteroskedastic and non-gaussian. This is confirmed by model fit and diagnostics measures. Models M1d and M1e seem to perform the best. Since the differences between the predictions from these models are minimal, we will choose the simpler model, i.e. M1d. Performance of various models measured in the training set (1a) and the test set (1b); According to a log. Score and normality of pseudo-z-scores. Avg. log. Score: bigger is better, skewness: closer to 0 is better (no skewness), kurtosis: closer to 0 is better, higher than 0 implies heavy tails. W: is a test statistic from the Shapiro-Wilk normality test. Closer to 1 means data pseudo-z-scores are more normal. Both AIC and BIC chose the best model.

M1a M1b M1c M1d M1e M1a M1b M1c M1d M1e -5422.29 0 -5429.8 -5444.86 0

0 -5857.19 -5858.21 (detrended Q-Q) plot. The closer the points are to the horizontal line, the better calibrated the model is. The confidence band shows where 95% of the points from a perfectly calibrated model are expected to lie. We see a small non-linear and heteroskedastic effect in the data. According to the wormplot in panel B, models m1d and m1e have a better fit than alternative models but are very similar, seen as mostly overlapping points.

Here, we will examine different ways how the effect of sex on the cortical thickness can be modeled. We can choose not model sex at all (M2a, same model as M1d, but fitted in a sample containing both men and women), choose to model separate intercepts for men and women, while keeping the trajectories the same (m2b), choose to model different age trajectories for men and women (m2c), and last, choose to model all aspects of the model separately for men and women, which is the same as having two separate models (m2d): M2a: ℎ ∼ (µ, σ, ν, τ), µ = M2b: ℎ ∼ (µ, σ, ν, τ), µ = M2c: ℎ ∼ (µ, σ, ν, τ), µ = σ = (), ν = β σ M2d: ℎ ∼ (µ, σ, ν, τ), µ = σ =

( * ) + σ,1 () + β (), ν = β , τ = β µ,1 We could also make more complicated models by, for example, modeling the effects of age on ν and τ parameters, but we showed in the previous section that this is not beneficial in this case. Results can be seen in table 2. The results show that sex effects are minimal. In the next section, we are therefore not going to model specific sex effects.

All models performed similarly out of sample, indicating a very small sex effect in this data modality.

To showcase modeling of site-effects, we use a combined multisite dataset as in (Kia et al. 2020) . This dataset includes data from various publicly available studies. The processing details can be found in Kia et al. (2020) . As in previous examples, we will model median cortical thickness data obtained using freesurfer 5.3 and 6.0. For simplicity, and due to results from previous section, we do not explicitly model sex effects, although combining modeling of site-effects and sex effects is also possible.

We compare different ways of modeling effects of scan-site. We can choose not to do anything (M3a), model site intercepts for µ as a fixed effect (m3b) or a random effect (m3c), and model site effect for µ and σ as fixed (m3d) or random effects (m3e). We could also model site-specific smooth terms, which could be potentially pooled towards a common trend as in random effect modeling of intercepts see (Pedersen et al. 2019) . M3a: ℎ ∼ (µ, σ, ν, τ), µ = µ µ σ ν, τ = β

τ () + β, σ = σ (), ν = β ν, τ = β τ Results can be seen in table 3. We can see that modeling site effects for µ and σ both significantly improve predictions in the new subjects. This can be seen in the deviance, but also when looking at z-scores stratified by site. We can also see that if the site effects are modeled, there is a smaller site effect present in the pseudo-z-scores in the test set. There will always be some site effect noticeable in the test set unless the sample size is huge, because, just by chance, some site effects will be overestimated or underestimated. Also, due to the sampling variability in the test set, the differences in averages and standard deviations of pseudo-z-scores across sites cannot be exactly zero. model logscore skewness kurtosis W