Correlation Constraints for Regression Models: Controlling Bias in Brain Age Prediction

<jats:p>In neuroimaging, the difference between chronological age and predicted brain age, also known as<jats:italic>brain age delta</jats:italic>, has been proposed as a pathology marker linked to a range of phenotypes. Brain age delta is estimated using regression, which involves a frequently observed bias due to a negative correlation between chronological age and brain age delta. In brain age prediction models, this correlation can manifest as an overprediction of the age of young brains and an underprediction for elderly ones. We show that this bias can be controlled for by adding correlation constraints to the model training procedure. We develop an analytical solution to this constrained optimization problem for Linear, Ridge, and Kernel Ridge regression. The solution is optimal in the least-squares sense i.e., there is no other model that satisfies the correlation constraints and has a better fit. Analyses on the PAC2019 competition data demonstrate that this approach produces optimal unbiased predictive models with a number of advantages over existing approaches. Finally, we introduce regression toolboxes for Python and MATLAB that implement our algorithm.</jats:p>

section_content: As the world's population ages, early detection and prevention of neurological aspects of aging, such as cognitive decline and dementia, is a public health priority and challenge. Pathological aging could be indicated by the level of deviation from the typical pattern of aging in healthy individuals (1). There has been growing interest in developing statistical approaches in order to identify individuals deviating from a healthy brain aging trajectory (2). To this end, a metric referred to as brain age delta, defined as the difference between brain-predicted age and chronological age, has been proposed as an index of the level of neuropathology in aging (2) (3) (4). Investigating the association between this metric with demographics, and lifestyle and cognitive variables can deepen the understanding of the processes that underpin healthy aging (5). In clinical research, brain age delta has the potential to index the severity of premature aging in patients suffering from disease. Among others, a higher delta has been associated with lower fluid intelligence and higher mortality (1), risk for developing Alzheimer's disease (6), severity of schizophrenia and depression (7). 

Establishing a good estimate of brain age delta is faced with important methodological challenges. The first challenge relates to the kinds of features (biomarkers) and predictive models that are used to build a brain-age model. A number of brain metrics have been considered as features for regression models; for example, structural networks (8), cortical thickness (9), functional connectivity patterns (10, 11), and raw T1-weighted images (1, 2). Likewise, a variety of regression models has been tested, from linear regression models such as lasso and support vector regression (10, 11) to convolutional neural networks [CNNs; (2, 8) ]. The quest for more accurate brain-age models lies at the heart of the Predictive Analytics Competition (PAC) 2019 upon which the eponymous Frontiers Research Topic is founded 1. 

A more fundamental methodological challenge, and the starting point for this paper, is the very operationalization of the brain age delta. If we denote the chronological age for a set of participants as a vector y, the ages predicted on the brain scans as ŷ, and the residuals as e = y -ŷ, then the negative residual δ = -e (i. e., predicted brain age minus chronological age) is usually defined as the brain age delta. This metric has been shown to be problematic. A predictive bias manifesting as an overprediction of the age of young individuals and an underprediction for elderly individuals has led to much speculation and investigation (2, 3, 9, (12) (13) (14) (15). 

A useful quantification of this effect is the correlation between chronological age and delta, corr(y, δ), which we will refer to as age delta correlation (ADC) in the rest of the paper (see Figures 1A, B ). An analysis by (15) showed that negative ADC is ubiquitous across a range of aging datasets and regression models and independent of the age range included in the dataset. A theoretical analysis by (14) showed that this effect is an inevitable property of regression, further aggravated by regression dilution (16, 17), and hence not limited to aging data. The potential danger of non-zero ADC lies in spurious associations with other covariates: brain age delta can be trivially correlated with demographic or cognitive variables if the latter are correlated with chronological age as well. Le et al. (14) found that associations between residuals and variables obtained from clinical interviews and neuropsychological testing largely disappear when residuals are corrected for chronological age. To avoid these spurious correlations, several authors have suggested following up the regression analysis with a correction step wherein the effect of age is removed from the residuals (1, 3, 12, 14). Brain age delta is then calculated by the following two-stage approach: Despite the significant methodological progress that has been made there are still concerns that warrant attention. First, the correction approach is an ad hoc fix because the models' predictions do not take ADC into account. Second, in a strictly sequential two-stage approach it is not clear whether the resultant brain age delta is optimal in terms of predictive accuracy. Both issues can be addressed when prediction and correction are unified within a model. Third, in predictive settings with training and test set, zeroing ADC on the training set is not of primary importance. Rather, it would be useful to have a model that is able to finely control the trade-off between ADC and predictive accuracy in order to optimize its performance on the test set. 

The aim of this study was to address these points by introducing modifications to three different regression models (Linear, Ridge and Kernel Ridge regression). Our models explicitly control for ADC on the training set without the need for an additional correction after model training. This was realized by formulating constrained optimization problems that incorporate age delta correlation as additional constraints. Our approach offers the following features: We sampled 100 random Gaussian and sigmoid datasets as specified before. We created random regression models by sampling slopes and intercepts from a uniform distribution (the sampling space included the OLS coefficients). When calculating ADCs on random models, a clear bias toward large negative ADC values is evident. Its distribution is more peaked for the Gaussian data. (E) ADC density plot. Simulated data (1,000 samples, 5 features, linear function, Gaussian noise) is split into 50% train, 50% test data. An OLS model is calculated on the train data and ADC is calculated on both train and test data. For the test data, the variability of the ADC is larger than for train data, and its mode is slightly more negative. 

the entire operation of the model is represented by its regression coefficients and quantities derived from them (18) (19) (20). Furthermore, we show in Section 2. 7 that these quantities are not affected by the choice of the correlation bound (a hyperparameter in our model). 

Some of the existing approaches share some of the listed features. For instance, whereas in (14) test set residuals are used for correction, Beheshti et al. (12) learns the correction parameters from the training set and applies them to the test set, in line with good practice for predictive models. However, to the best of our knowledge, this is the first study to prove optimality of the models and introduce "soft" correlation bounds for fine control of ADC.
section_name: INTRODUCTION

section_content: The section 2 is organized as follows. In Section 2. 2, we introduce Linear regression, Ridge regression and Kernel Ridge regression. 

In Section 2. 3, we review existing ways to quantify brain age delta. In Section 2. 4, we revisit the mathematical basis for ADC in the context of the three regression models. In Section 2. 5 we develop our approach by adding correlation constraints to the model training stage that allow for a precise control of ADC. In Section 2. 6 it is shown that the brain age estimates obtained with our models are closely related to existing correction approaches. In Section 2. 7, we investigate how to interpret the models. In Section 2. 8, we introduce corresponding toolboxes for Python and MATLAB that implement them. Finally, in sections 2. 9 and 2. 10, we describe our analysis of the PAC2019 competition data.
section_name: METHOD

section_content: Table 1 defines the most important mathematical symbols used in the paper. Whenever the symbol represents a vector or matrix, its dimensionality is given in the second column. In general, matrices are denoted as uppercase boldface symbols (e. g., X, H), vectors as lowercase boldface symbols (e. g., β, x), and scalars as lowercase normal face symbols (e. g., n, ȳ). 

Note that we use terminology common in the machine learning literature. In particular, features are also known as predictors or independent variables in the regression literature, the vector of targets (chronological age) is also known as response vector or dependent variable, and training is also known as fitting. For standard regression models, regression coefficients are denoted as β. If an intercept is included in the model, we assume that a column of 1's is added to the matrix of features X. Sometimes we explicitly denote the intercept as β 0 and the non-intercept coefficients as β 1 : p. For models with correlation constraints, we use the notation b for the regression coefficients with b 0 and b 1 : p defined analogously.
section_name: Notation

section_name: Regression Models

section_content: Ordinary least-squares regression, often just called linear regression, can be formulated as a set of n equations of the form

where y i is the i-th response, x ij is the jth predictor value in the i-th sample, the β's are regression coefficients, and ǫ ∼ N (0, σ 2 ) is an error term. Using matrix notation, this set of equations can be written more succinctly as:

where y = [y 1, y 2,. . ., y n ] ⊤ comprises the responses, X ∈ R n×(p+1) is the matrix of features including a column of 1's for the intercept term, β ∈ R p+1 is the vector of regression coefficients and ǫ = [ǫ 1, ǫ 2,. . ., ǫ n ] ⊤ collects all error terms. Training a model implies finding an estimate for β such that y ≈ Xβ. In OLS regression this is achieved by minimizing the sum of squared errors ||y -Xβ|| 2 2. Denoting ŷ : = Xβ this can be formulated as the unconstrained optimization problem minimize ||y -ŷ|| 2 2. 

(

The solution is given by
section_name: Ordinary Least-Squares (OLS) Regression

section_content: OLS regression suffers from high variance and does not have a unique solution if the number of the features p is larger than the number of samples n (21, 22). Ridge regression is a regularized version of OLS regression that is useful for data that suffers from multicollinearity. The model is regularized by adding an ℓ 2 penalty that shrinks the weights toward zero. For a given regularization parameter λ ≥ 0, Ridge regression can be formulated as the unconstrained optimization problem

The first term is the least-squares term from Equation ( 2 ). The second term penalizes elements of β from becoming too large. For λ = 0 Ridge regression reduces to OLS regression. The solution is given by

where p+1 ) is an identity matrix. Since the intercept term is not regularized, I is modified such that the 1 in the first row/column is replaced by a 0.
section_name: Ridge Regression

section_content: Other variants of linear regression e. g., lasso (23) and elastic net (24), do not have a closed form solution but rely on iterative optimization, so they do not lend themselves to the analytical approach developed in this paper.
section_name: Other Linear Regression Models

section_content: A non-linear version of Ridge regression can be developed by applying a non-linear transformation to the features and then performing Ridge regression on these transformed features (25). Let this transformation be represented by a map φ : R p → F from input space to a higher-dimensional Reproducing Kernel Hilbert Space and (X) = [φ(x 1 ), φ(x 2 ),. . ., φ(x n )] ⊤ (26, 27). The solution is given by replacing X by (X) in Equation ( 3 ),

This solution, also known as primal solution, is of limited practical use, since the feature space is often too highdimensional to represent β krr and (X). As an alternative, the convex optimization problem can be rewritten into its dual Lagrangian form first (28). The resultant dual solution is given by

The equivalence between the primal and dual solution can be verified by left-multiplying both solutions with ( (X) (X) ⊤ + λI). Since β krr cannot be represented directly, we instead calculate a vector of dual weights α ∈ R n. To this end, define K = (X) (X) ⊤ as the kernel matrix with K ij = k(x i, x j ) for a kernel function k. Then the vector of dual weights is given by the latter part of Equation ( 4 ), α = (K + λI) -1 y. 

(

Using α the predicted response to a test sample x can be rewritten in terms of kernel evaluations:
section_name: Kernel Ridge Regression (KRR)

section_content: The regression of age on brain features often leads to a biased model that manifests as an overprediction of the age of younger individuals and an underprediction of the age of elderly ones. This effect can be quantified as a negative age delta correlation (ADC) denoted as corr(y, δ). In the literature, ADC has been set to zero by adding a second stage to the analysis wherein the regression predictions from the first stage are corrected. Hence, brain-age delta prediction can be formulated as the following two-stage approach:

(a) Brain Age Prediction. Train a regression model f to predict age such that y ≈ f (X). The negative residuals, denoted as

represent the uncorrected brain age delta. (b) Correction of Brain Age Delta. A number of authors proposed correction procedures to rid δ of ADC (1, 3, 4, 12, 14). Many of these approaches are mathematically equivalent. They boil down to two approaches that yield two differently corrected residuals δ 1 (approach 1) and δ 2 (approach 2). These two approaches are discussed in detail in the following two subsections.
section_name: Calculating the Brain-Age Delta

section_content: Approach 1 has been proposed by a number of authors (3, 12, 14) and boils down to the following operation: Train a simple regression model δ ≈ yβ 1 + β 0 to remove the linear effect of age from the delta estimate. The new estimate

represents corrected brain age delta which is uncorrelated with age. We can inspect the value of β 1 by taking the simple linear regression formula: β 1 = r δy s δ s y, where r δy is age delta correlation (ADC) which is negative (Section 2. 4). s δ is the standard deviation of the residuals and unnormalized square root of the residual sum of squares (RSS), s y is the standard deviation of the responses and the unnormalized square root of the total sum of squares (TSS). In OLS, Ridge, and Kernel Ridge regression, we have TSS ≥ RSS. Hence, β 1 ∈ [-1, 0]. Combining Equations ( 8 ) and (7) the entire model can be written in one equation as

where (1 + β 1 ) ∈ [0, 1]. From Equation ( 9 ), we can see that the correction does not affect the predictions f (X). Instead, it implies shrinking y. 

This approach is invalid in a predictive modeling framework because it corrects the data y, not the predictions f (X). Beheshti et al. (12) report a lower error and a larger R 2 value compared to approach 2 introduced in the next section. However, since this effect is obtained by shrinking the data it can be considered as an artifact of this approach.
section_name: Approach 1: Scale Down y (Chronological age)

section_content: As an alternative approach, it has been suggested that ŷ instead of δ should be used in the regression (1, 4). To this end, train a simple regression model ŷ ≈ yβ 1 + β 0. Then define corrected predictions ŷ2 as

with corresponding brain age delta δ 2 = ŷ2y. This modified age delta estimate is again uncorrelated with age. As in approach 1, the correction is performed using simple linear regression and we have β 1 = r ŷy s ŷ s y. s ŷ is the standard deviation of the predictions and the unnormalized square root of the explained sum of squares (ESS). In OLS, Ridge, and Kernel Ridge regression, we have TSS ≥ ESS and r ŷy ∈ [0, 1]. This implies that β 1 ∈ [0, 1]. Combining Equations ( 10 ) and (7) this can be written in one equation as

Comparing Equations ( 9 ) and ( 11 ), we see that in approach 1 the data vector y is scaled down whereas in approach 2 the predictions f (X) are scaled up. In Section 2. 6, we show that the two types of corrected residuals are actually identical up to scaling and therefore corr(δ 1, δ 2 ) = 1. Consequently, they perform equally well on secondary analyses e. g., relating brain age delta to cognition. They are further closely related to our zero correlation constraint (Section 2. 5. 1). In a predictive modeling framework, we consider approach 2 as preferable since corrections should be applied to the model not to the data.
section_name: Approach 2: Scale up f(X) (Predicted age)

section_content: The theoretical basis for negative ADC has already been discussed in (14). In particular, the authors highlighted that ADC ≤ 0 for any sensible regression model. Here, we discuss ADC more specifically for the three regression models introduced above. We start with OLS regression. Let us expand the age delta correlation term as corr(y, δ) = y ⊤ (ŷy) ||y|| ||ŷ -y|| (12) where to simplify the notation we assume that X and y have been centered. The sign of the correlation is determined by the numerator. Defining H = X(X ⊤ X) -1 X ⊤ and writing ŷ = Hy we obtain

where the inequality follows from the fact that I -H is symmetric and idempotent and therefore positive semi-definite (29). This implies that the ADC is always non-positive in OLS regression. Ridge regression coincides with OLS regression for λ = 0. As λ increases, β tends to zero due to the shrinkage effect of the regularization (21). This implies that ŷ → 0 and therefore corr(y, δ) → -1 as λ increases. This is illustrated empirically in Figure 1C. The same argumentation holds for Kernel Ridge, one only has to replace X by (X). Often, Kernel Ridge models will have a smaller prediction bias because their higher complexity allows for a better fit to the data. Furthermore, Figure 1D shows a clear bias toward large negative ADC values when regression coefficients are randomly sampled. Figure 1E shows that the bias persists in both train and test sets. Together, these results suggest that a negative ADC is inevitable and that regularization further exacerbates this effect, in line with previous work (3, 14).
section_name: Negative Age Delta Correlation (ADC)

section_content: The regression problems defined in Section 2. 2 can be cast as unconstrained optimization problems. The optimization involves the minimization of a loss function L which measures the amount of discrepancy between the true responses y and the model predictions ŷ:

In OLS regression the loss function is the squared loss

2 in Ridge and Kernel Ridge regression. To control for age delta correlation in the training data, we can add a correlation constraint that caps the permitted magnitude of correlation between the brain age delta and age. To this end, consider the constrained optimization problem minimize L(y, ŷ)

The same loss function as before is minimized. However, the set of feasible solutions is limited to solutions for which the absolute value of the correlation does not exceed ρ, where ρ ≥ 0 is the correlation bound selected by the user. As a special case of Equation ( 15 ), we can consider the case ρ = 0, that is, the responses have to be perfectly uncorrelated with the residuals: minimize L(y, ŷ) subject to corr(y, δ) = 0. 

(16)

We will address the latter case first and see that it leads to a simple solution. In the following, we will assume that X and y have been centered and the model contains no intercept since this simplifies the equations. This does not limit the generality of our results. As shown in the Supplementary Material (Section 1), a model with centered data and without intercept yields the same regression coefficients as a model with intercept. In other words, we can first calculate β 1 : p on the centered data and subsequently calculate the intercept β 0 to obtain the model for non-centered data.
section_name: Correlation Constraints for Regression

section_content: A hard correlation constraint can be set that requires the correlation between the residuals and the response values to be zero, that is |corr(y, δ)| = 0. In the Supplementary Material (Section 2) the optimal solution, b, is derived for OLS, Ridge, and Kernel Ridge regression. It is given as a scaled version of the standard, unconstrained solution

where β is the standard OLS, Ridge, or Kernel Ridge solution and it is assumed that X and y have been centered. Using Equation ( 6 ), we can see that for Kernel Ridge regression this translates into a scaling of the dual weights

The scaling factor θ 0 is given by

with model-specific hat matrices H:

Intercept term: If the data has not been centered and the model includes an intercept term, y and y ⊤ Hy need to be centered before calculating θ. The intercept b 0 can be obtained from the equation
section_name: Zero Correlation Constraint

section_content: The zero correlation solution successfully removes the correlation between residuals and responses. However, it does so at the cost of goodness of fit to the training data. Furthermore, in predictive modeling, perfect control of ADC on the training set is less important than good predictive performance and low bias on the test set. This suggests the need for a softer constrained optimization solution wherein the equality constraint is replaced by an inequality constraint. In the Supplementary Material (Section 3) it is shown that the optimal solution is again given by scaling, b 1 : p = θ ρ β 1 : p, where there are now two possible solutions for the scaling factor,

where c = (y ⊤ ŷ) 2ρ 2 ||y|| 2 ||ŷ|| 2 and ŷ is the predictions under the unconstrained OLS, Ridge, or Kernel Ridge model. The two solutions define an interval [θ 1, θ 2 ]. Setting θ to any value within this interval will guarantee -ρ ≤ corr(y, δ cc ) ≤ ρ, where δ cc is the brain age delta under the correlation constrained models. Setting θ = θ 1 or θ = θ 2 will set corr(y, δ cc ) = -ρ or corr(y, δ cc ) = ρ, respectively. From this we can derive the following algorithm for correlation constrained models with an inequality constraint:

1. Calculate the standard, unconstrained solution for the model (OLS, Ridge, or Kernel Ridge). If |corr(y, δ)| ≤ ρ, the unconstrained solution does not violate the correlation constraints so we are done. 2. If |corr(y, δ)| > ρ, calculate θ 1,2 and set θ to the value that smaller in absolute value. This will assure that corr(y, δ cc ) = -ρ if corr(y, δ) < -ρ. 

Figure 2 depicts the geometrical intuition underlying the correlation constraints. Without constraints, the solution is the minimum of a quadratic function (Figure 2A ). For a zero correlation constraint, the set of feasible solutions is reduced to a line within this space (Figure 2B ). For a bounded correlation constraint, the set of feasible solutions is the space between two paraboloids (Figure 2C ). In both cases, the correlation constraints lead to a larger slope for the regression coefficients as compared to the unconstrained model (Figure 2D ).
section_name: Bounded Correlation Constraint

section_content: In Section 2. 3, we reviewed the two main approaches for correcting brain age delta used in the literature. Here, we investigate their mutual relationship as well as their relationship to our approach. Without loss of generality we assume that y and X have been centered (see Supplementary Material, Section 1). Let us start from Equation (8) corresponding to approach 1. The regression slope β 1 for a simple linear regression model is given by

Writing δ = ŷy and expanding the correlation term as in Equation (S5) (Supplementary Material), we find that

with θ 0 as defined in Equation (19). Therefore, the solution to Equation ( 8 ) is given by

Alternatively, using approach 2 (correction of predictions) and starting from Equation (10) we perform a regression of ŷ on y. Again, this is a simple linear regression model whose slope is given by

yielding the corrected predictions ŷ2 = ŷ/β 1 = θ 0 ŷ with corresponding brain age delta

This solution uses a scaling of ŷ by θ 0 and is thus equivalent to our zero correlation solution when an OLS model is used. Furthermore, comparing Equations ( 20 ) and (21) we see that both solutions are proportional to each other. Their relationship is given by δ 2 = θ 0 δ 1 (22) and therefore corr(δ 1, δ 2 ) = 1. In other words, the brain age delta estimates from the two approaches used in the literature are identical up to a scaling factor of θ 0.
section_name: Relationship Between Zero Correlation Constraint and Existing Correction Approaches

section_content: Maximizing predictive performance is the primary objective when optimizing statistical models. To empirical researchers, understanding what the regression model learns from the data is useful, too. Linear regression models such as OLS and Ridge can be interpreted in terms of their β coefficients. To keep the notation simple, we will assume that the data has been demeaned and the model contains no intercept term. If the data is standardized, large components of the β can be interpreted as features that are relevant to the regression task. Since our models combine prediction and correction into a single task, the β's capture the entire operation of the model. Importantly, the choice of the correlation bound ρ does not change the interpretation. Since the regression coefficients in our models are just scaled versions of the original regression coefficients, b = θ β for some θ ∈ R, the choice of ρ does not affect the ratio between any pair of coefficients. 

For collinear data, coefficients can become uninterpretable with large weights for features that are not related to the target variable. In this case, structure coefficients (18) and activation patterns (19) have been proposed as alternatives metrics. An activation pattern is given by

where a ρ ∈ R p is the activation pattern and is the data covariance matrix. Let a β = β be the activation pattern for the standard (uncorrected) model. Then setting ρ merely scales the activation pattern by θ ρ since a ρ = b = θ ρ β = θ ρ a β. An example for an activation pattern is depicted in Figure 5 and discussed in Section 3. 3. 

Structure coefficients are given by the vector of Pearson correlations between ŷ = Xb and each column of X, where the i-th column is denoted as X i. Since the correlation coefficient is invariant to constant shifts and scaling, we have corr(X i, Xb) = corr(X i, Xβ), that is, structure coefficients are invariant to the choice of ρ. 

For kernel methods, these coefficients are generally not available. As an alternative, we considered the vector of partial derivatives of regression model with the respect to each of the features,

where k is a kernel function (20) and i refers to the index of the feature. Using Equation ( 18 ), we see that for a correlation constrained Kernel Ridge model f ρ these partial derivatives are again just scaled versions of the unconstrained version

To summarize, our models are interpretable because the regression coefficients capture the whole operation of the model. Furthermore, the choice of the correlation bound does not affect the interpretation of the model: structure coefficients are invariant to the choice of ρ, whereas regression coefficients, activation patterns and partial derivatives are merely scaled by a constant factor.
section_name: Interpretability

section_content: The Linear, Ridge, and Kernel Ridge regression models with correlation constraints presented in this paper have been implemented in Python and MATLAB. For Python, the models are available on GitHub (github. com/treder/correlationconstrained-regression). They inherit from and are fully compatible with the Scikit-Learn framework (30). The models extend Scikit-Learn's LinearRegression, Ridge, and KernelRidge models with an additional parameter correlation_bound that corresponds to ρ in Equations ( 15 ) and ( 16 ). Setting the parameter to 0 enforces a zero correlation constraint whereas setting it to a positive value bounds the correlation accordingly. For MATLAB, the models have been integrated into MVPA-Light (31), an open-source machine learning toolbox. By setting the hyperparameter correlation_bound, ADC can be controlled in the same way as for the Python-based models. Code examples for both Python and MATLAB can be found on the GitHub page.
section_name: Toolbox

section_content: Neuroimaging data supplied within the Predictive Analytics Competition were fully pre-processed T1-weighted MRI scans from 2,640 training set and 660 validation set subjects as described previously (2). All normalized 3D maps of gray matter (GM) and white matter (WM) volume were used to create group GM and WM masks. Each GM and WM image was smoothed using an 8-mm Gaussian kernel, masked and concatenated into a vector of 153,237 and 86,143 voxels, respectively. 

Concatenated GM and WM images were intensity normalized and submitted to Independent Component Analysis using the Group ICA of fMRI toolbox [https://trendscenter. org/software/ gift/; (32) ]. The optimal number of components of the ICA decomposition (72 and 99 for GM and WM images, respectively) was determined using Principal Component Analysis (PCA) with minimum description length (MDL) model order selection criteria (33). Normalized features in the training model, based on 2,640 participants, included scores for all GM and WM components (N = 171). Additional covariates included total GM, total WM, gender and dummy coding for 17 scanning sites. 

We also included low-order interaction terms (such as bivariate interaction) between total GM, total WM, gender, PC1 and PC2 scores.
section_name: Neuroimaging Data

section_content: We performed brain age prediction using Python with Scikit-Learn and our custom extensions. The models were tested in three different conditions: Unconstrained (using Scikit-Learn's models without correlation constraints), zero correlation (ρ = 0) (using our extensions with a correlation constraint of 0), bounded correlation (ρ = 0. 1, 0. 2, 0. 3) (using our extensions with a correlation bound of 0. 1, 0. 2, 0. 3), and approaches 1 and 2 from the literature introduced in Section 2. 3. To obtain both insample and out-of-sample statistics, models were applied to both training and test data. To estimate the variability of predictive performance, we performed 100 iterations. In every iteration, training data was randomly sampled from the training set using bootstrapping. The 171 Independent Components were used as features. 

Three regression models were considered, OLS, Ridge regression, and Kernel Ridge regression with a RBF kernel. For both Ridge and Kernel Ridge regression, hyperparameters were tuned using a grid search with Scikit-Learn's GridSearchCV and five-fold cross-validation. For Ridge regression, the regularization parameter was tuned using candidate values α = (10 -3, 10 - 2, 10 -1, 1, 10). For Kernel Ridge with a RBF kernel, kernel width γ = (100, 10, 1, 10 -1 ) and α = (10 -3, 10 -2, 10 -1, 1, 10) were tuned. The resultant best model was used to calculate in-and out-of-sample metrics. Mean absolute error (MAE) and age delta correlation (ADC) served as metrics. Denoting train and test sets as T R and T E, MAE was estimated as

where f is a regression model trained on the training data and |T R| and |T E| are the sizes of the train and test sets, respectively. Similarly, ADC was calculated separately for the predictions in the train and test sets. All analyses were performed on a Desktop computer with an Intel Core i7-6700 @ 3. 40 GHz x 8 CPU with 64 GB RAM running on Ubuntu 18. 04. The analysis code is available as part of the toolbox 2. (standard regression model), bounded correlation constraint (ρ = 0. 1, 0. 2, 0. 3), zero correlation constraint (ρ = 0), and approaches 1 and 2 from the literature introduced in Section 2. 3. The same data is presented in tabular form in Table 2. Bonferroni correction was used in case of multiple comparisons.
section_name: Brain Age Prediction

section_name: RESULTS

section_content: With respect to mean absolute error (MAE), a significantly lower error was found on the train data compared to test data (Wilcoxon signed-rank test, W = 62, 449, p < 0. 0001). For the constrained models, we performed a regression analysis of MAE on ρ. We found significant negative slopes for all models and train and test sets separately (all p < 0. 0001), indicating that MAE decreases significantly as ρ increases. On the test data, we used Wilcoxon signed-rank tests to compare MAE for the unconstrained model with MAE for each of the constrained models. A significantly higher MAE was obtained for most unconstrained models (all p < 0. 0001). The only exception was approach 1 ("scale down y") wherein the relationship is reversed: MAE decreases after correction. This is an artifact of the fact that the correction is applied to the data, not the predictions (see Section 2. 3). 

With respect to age delta correlation (ADC), our models perfectly controlled for ADC on the training data. ADC was significantly larger in magnitude for the test data compared to the train data (W = 877, 894, p < 0. 0001). Linear regression of ADC on ρ showed a significant negative slope for all models on both train and test sets (all p < 0. 001). On the test set, ADC was significantly larger in magnitude for the unconstrained model than for the constrained models (all p < 0. 0001), suggesting that the constrained models also control ADC on the test set.
section_name: Brain Age Prediction

section_content: To better characterize how the choice of ρ mediates the tradeoff between ADC and MAE, we repeated the prediction analysis. This time the correlation bound was varied in small steps of 0. 02. Results averaged across 100 bootstrap iterations are depicted in Figure 4. In line with the results above, MAE generally decreases with increasing ρ and eventually flattens off. ADC decreases roughly linearly with ρ. It flattens off at a value of ρ that corresponds to the ADC value of the uncorrected model, since no correction needs to be applied when ρ > |ADC|. Furthermore, MAE and ADC change similarly for train and test set, albeit with different slopes. This is a useful observation since it suggests that the hyperparameter ρ can be optimized on the training set alone, in line with good practice in predictive modeling.
section_name: ADC-MAE Trade-Off

section_content: Figure 5 shows an activation pattern for an OLS model with ρ = 0 trained on the whole training set, using gray matter Independent Components (ICs) only. To create the maps, a vector of activation patterns a was calculated for the gray matter ICs (see Section 3. 3). Since each IC corresponds to a brain map, we multiplied each entry of a with its corresponding map and added up the maps. The figure depicts the resultant summed map indicating an age-related decrease in intensity values in deep cortical areas. An age-related increase in intensity values was observed at the boundaries between gray matter and other tissue types, likely reflecting CSF signals in older adults (34). Crucially, the map is independent of the choice of ρ, only its scaling is affected by ρ. This illustrates that the interpretation of the models is not affected by a change of ρ.
section_name: Interpretability via Activation Patterns

section_content: A predictive bias manifesting as an overprediction of the age of young individuals and an underprediction of the age of elderly individuals has been consistently reported in the brain age literature (2, 3, 14, 15). It can be quantified as age delta correlation (ADC), that is, the correlation between brain age delta (predicted age minus chronological age) and chronological age. We introduced modifications to three popular regression models, OLS, Ridge and Kernel Ridge regression, that effectively control ADC. To this end, we introduced a hyperparameter ρ that can be set by the user. It represents a correlation bound that controls the maximum permissible ADC. The resultant models are optimal in that they give the lowest mean-squared error on the training set (among all models from the same class) while controlling for ADC. Our models were tested on the PAC competition data using different values for ρ. The models not only perfectly controlled ADC on the training data, they also approximately controlled ADC on unseen test data (see Figure 3 ). For all constrained models, ADC on the test set was lower than for the unconstrained models. The flip side of this was an increase of mean absolute error (MAE) for our constrained models as compared to the unconstrained model, but often this increase was modest. For instance, for an OLS model MAE increased from 5. 31 for the unconstrained model to 5. 35 for the model with ρ = 0. 3, an increase of only 0. 8%. Across all models in the test set, we found that an increase in ρ led to a decrease in MAE. This suggests that ρ can be used as a lever to finely control the trade-off between predictive performance (MAE) and age delta correlation (ADC). 

In the same analysis, we included the two existing correction methods used in the literature, denoted as approach 1 (3, 12, 14) and approach 2 (1, 4) discussed in detail in Section 2. 3. For the special case of using a OLS model with a zero correlation constraint, ρ = 0, our models' brain age deltas are equivalent to approach 2. In Section 2. 6, we furthermore show that the brain age deltas in approaches 1 and 2 are actually identical up to scaling. They differ only by the scaling factor θ 0 defined in Equation (19). In particular, in approach 1 chronological age is scaled down by this factor before calculating brain age delta, whereas in approach 2 f (X) (predicted age) is scaled up by the same amount. The downscaling in approach 1 leads to a lower MAE which is even smaller than for an uncorrected model. We would like to stress that this is because approach 1 corrects the data, not the predictions. This is not permissible in predictive modeling and approach 2 should be preferred. 

The hyperparameter ρ controlling ADC has to be selected by the user. Figure 4 shows how the choice of ρ affects both MAE and ADC. A possible selection criterion would involve defining a maximum permissible MAE or ADC value and choosing ρ accordingly. A more comprehensive analysis would take into account any follow-up analyses. For instance, brain age delta is often correlated with cognitive variables in a second step. An optimal selection of ρ could involve e. g., the regression slope or pvalue of the association on the training set. A detailed exploration of how the choice of ρ affects follow-up analyses requires a dataset including cognitive test scores and is left for future work. 

In terms of interpretability, our models offer a greater degree of transparency than the traditional two approaches because the model coefficients capture the entire brain age prediction pipeline (i. e., both prediction and correction). Moreover, in Section 2. 7 we show that the interpretation is not affected by changes in ρ, a hyperparameter in our models. Structure coefficients (18) are invariant to changes in ρ, whereas the other metrics are simply scaled by a constant value. 

A limitation of our study is that it only covers OLS, Ridge, and Kernel Ridge regression. In the age of deep neural networks, the focus on linear and kernel methods may seem very limiting. Therefore, we would like to emphasize that linear regression models lie at the heart of many non-linear approaches including Convolutional Neural Networks. As illustrated in Figure 6, nonlinear approaches can often be conceived of as linear regression operating on non-linearly extracted features. Linear models can uniquely combine predictive power with computational efficiency and interpretability. In line with this, (35) found that kernel regression was as performant as deep neural networks when predicting phenotypes from functional connectivity data. Nevertheless, future work could address the incorporation of correlation constraints into other models classes such as Lasso (36), Support Vector Regression (37) or CNNs. Since these models use iterative optimization, a possible approach could be adding the correlation term -corr 2 (y, δ) directly to the loss function. Alternatively, since y is constant this can be simplified to the quantity y ⊤ ŷ/||δ|| 2. 

On a more speculative note, future development of the brain age delta metric might benefit from work on errors-invariables models (38) (39) (40) or measurement error models (16). Standard linear regression models assume that chronological age has been measured with an error whereas the brain data is noisefree. It is more likely that the opposite is true: chronological age can be measured with high accuracy but there is noise and individual variability in the brain scans. Not accounting for measurement error in the features leads to regression dilution which in OLS regression manifests as an underestimation of the regression coefficients. This phenomenon is known in the brain age literature (3, 14). Our scaling factor θ inflates the regression coefficients and therefore un-dilutes the model, but it is not clear to the authors whether there is a more formal relationship between correction of the residuals and measurement error. Unfortunately, estimating measurement error in brain scans requires repeated sampling which is often not available. 

Concluding, without accurate control for ADC, the use of brain age delta can lead to false associations with other phenotypes and limit our understanding of the processes that underpin brain aging. We highlighted the importance of estimating brain age delta and controlling for age delta correlation within a given model, as we introduced a novel class of regression models that allow for fine control of ADC. Our solution is optimal on the training set and shows approximate control of ADC on the test set. In an era of "big data" predictive modeling, this approach nicely dovetails with strategies to develop reliable models that generalize to independent test sets for use in personalized and precision medicine (41).
section_name: DISCUSSION

section_name: ACKNOWLEDGMENTS

section_content: We would like to thank Franz Király for insightful comments on the mathematical proofs. All authors contributed to the article and approved the submitted version.
section_name: ACKNOWLEDGMENTS

section_name: FUNDING

section_content: DS was funded by the SAMRC. KT was supported by the British Academy Postdoctoral Fellowship ( PF160048 ) and the Guarantors of Brain ( 101149 ).
section_name: FUNDING

section_name: DATA AVAILABILITY STATEMENT

section_content: The neuroimaging dataset analyzed in this study has been provided by the PAC 2019 team. All the program code and functions can be found in the following repository: https://github. com/treder/correlation-constrained-regression.
section_name: DATA AVAILABILITY STATEMENT

section_name: AUTHOR CONTRIBUTIONS

section_content: MT conceptualized the approach and developed the toolboxes for Python and MATLAB. KT processed the MRI data. MT and KT performed the age analyses. MT and JS performed the mathematical analysis. MT and KT wrote the manuscript with additional contributions, revision notes and comments from all authors.
section_name: AUTHOR CONTRIBUTIONS

section_content: The Supplementary Material for this article can be found online at: https://www. frontiersin. org/articles/10. 3389/fpsyt. 2021. 615754/full#supplementary-material
section_name: SUPPLEMENTARY MATERIAL

section_content: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
section_name: Conflict of Interest: