---
abstract: |
    Statistical shape analysis is a powerful tool to assess organ morphologies and find shape changes associated to a particular disease.  However, imbalance in confounding factors, such as demographics might invalidate the analysis if not taken into consideration.  Despite the methodological advances in the field, providing new methods that are able to capture complex and regional shape differences, the relationship between non-imaging information and shape variability has been overlooked. 

    We present a linear statistical shape analysis framework that finds shape differences unassociated to a controlled set of confounding variables.  It includes two confounding correction methods: confounding deflation and adjustment.  We applied our framework to a cardiac magnetic resonance imaging dataset, consisting of the cardiac ventricles of 89 triathletes and 77 controls, to identify cardiac remodelling due to the practice of endurance exercise.  To test robustness to confounders, subsets of this dataset were generated by randomly removing controls with low body mass index, thus introducing imbalance. 

    The analysis of the whole dataset indicates an increase of ventricular volumes and myocardial mass in athletes, which is consistent with the clinical literature.  However, when confounders are not taken into consideration no increase of myocardial mass is found.  Using the downsampled datasets, we find that confounder adjustment methods are needed to find the real remodelling patterns in imbalanced datasets. 
author:
- Gabriel Bernardino
- 'Oualid [Benkarim]{}'
- 'María [Sanz-de la Garza]{}'
- 'Susanna [Prat-Gonzàlez]{}'
- 'Alvaro [Sepulveda-Martinez]{}'
- 'Fàtima [Crispi]{}'
- 'Marta [Sitges]{}'
- 'Constantine [Butakoff]{}'
- 'Mathieu [De Craene]{}'
- 'Bart [Bijnens]{}'
- 'Miguel A. [González Ballester]{}'
bibliography:
- 'libraryMendeley.bib'
title: 'Handling confounding variables in statistical shape analysis - application to cardiac remodelling'
---

Introduction
============

Analysing the shapes of parts of biological organs and organisms has been the object of extensive study for over a century [@Thompson1942].  This interest in shape is also present in medicine: several studies have focused in the relationship between organ morphology and illness.  For instance, cardiac shape remodels to improve cardiac pressure/volume output under abnormal working conditions, and it is used to assess the presence/evolution of illness [@Arts1994; @Grossman1975].  In a nutshell, pressure overload produces concentric remodelling (thickening of the myocardium without dilation of the ventricle) to maintain wall stresses low, and volume overload dilates the ventricle without a myocardial mass thickening.  This is an oversimplication, as a volume overload will also increase pressure, and the exact remodelling mechanisms and triggers at a cellular level are still under research and discussion.  The classical way of analysing shape in the clinical community consists in manually extracting hand-crafted features, and analysing these shape descriptors.  These measurements are usually standardised and defined in guidelines [@Lang2015], and usually refer to global characteristics of the shape that carry little regional information. 

Nowadays, it is possible to acquire 3D images in clinical routine.  Furthermore, advances in computing permit to automatically segment the images and to generate personalised 3D models of the organs [@Ballester2000; @Mitchell2002; @Ecabert2006; @Bernard2018].  is a set of techniques to represent both shapes and images and do the analysis directly with these objects, not being limited to only analyse previously defined measurements.  is used in the medical imaging field, in order to identify and represent shape variability of those organs [@Cerrolaza2015; @Blanc2012; @Rajamani2007; @Sierra2006].  This allows expressing and quantifying regional shape patterns in a robust and objective manner, instead of working with a small set of predefined measurements on the shapes (like volumes and diameters).  can be used together with statistical learning techniques to construct models that find regional differences in anatomy that are associated to pathologies [@Zhang2014; @Singh2014a; @Varano2017; @Sarvari2017], based on a control and a pathological population.  Roughly, the typical framework consists of first using to construct an atlas of all shapes in the population, then use or another technique to find a low-rank representation of the shape space, and finally use a classification algorithm on that space to train a model that predicts the control/pathologic status. 

Beyond pathological variability, the shape of an organ also exhibits variability due to other factors, like lifestyle, gender, ethnicity, or size.  The framework described above uses the implicit hypothesis that differences in shape are only due to the pathology.  In some cases the remodelling associated to the pathology is prominent and easily identifiable.  However, in others, like subclinical or early-stage studies, differences can often be very subtle and less pronounced than demographic-related variability. .  Even in cases where the pathological remodelling is significant, and the populations are similar in terms of demographics, demographic-related variability will add noise to the analysis.  A number of authors have explored the usage of non-imaging information.  For instance Singh et al.  proposed a procedure similar to partial correlation between shape and several clinical variables while correcting for confounding variables[@Singh2014a].  Zhang et al adjusted by demographics in their studies of cardiac remodelling in myocardial infarction [@Zhang2014a], and Zhang et al and Mauger at al explored the relationship between shape and classical clinical measurements [@Zhang2017; @Zhang2014a; @Mauger2019].  However, not all authors include corrections for confounders, and their effect in shape analysis studies has not been yet quantitatively tested. 

In this paper, we present a framework to find differences between control and pathological populations that outputs the most discriminating shape pattern, that can be visualised for interpretability.  We quantitatively and qualitatively show the effect an imbalance in the confounding variables has in the analysis, and propose techniques to reduce that effect.  The proposed model consists of the following steps: (1) the construction of an atlas of the personalised 3D meshes automatically generated from images; (2) identification and removal of shape variability due to confounding variables; (3) dimensionality reduction and classification. 

To illustrate the framework, we use a dataset of cardiac involving sedentary controls and triathlon athletes.  This dataset was collected to study the remodelling due to the extended practice of endurance sport, which produces a volume overload to the heart.  This volume overload triggers compensatory mechanisms to improve cardiac output and withstand the increased pressure during exercise.  The whole of this remodelling is called the *Athlete’s Heart* and involves substantial changes in function and geometry at both rest and during exercise [@DAndrea2015; @Schiros2013a].  Although the remodelling is not yet completely understood, researchers have established a strong relationship between cardiopulmonary performance during exercise and cardiac geometry at rest [@LaGerche2012; @Scharhag2002]. 

Methodology
===========

The full process to compute the confounder invariant most discriminating shape pattern between two populations is summarised in Figure \[methodology:figure:framework\].  The presented framework consists of $3$ main steps:

1.   Compute the mean shape of the population, and register all shapes to this template. 

2.   Identify and remove shape variability attributable to confounding variables. 

3.   Train a classification model to obtain the most discriminating shape pattern between both populations and generate a visual representation of the most discriminating shape pattern

Atlas construction
------------------

From the sequence, we use a model-based automatic segmentation method to obtain personalised meshes with point-to-point correspondence [@Ecabert2006].  The method deforms a full heart (4 chambers) template mesh using a polyaffine deformation to match the myocardial boundaries.  The method includes slice correction to remove misalignment between consecutive slices. 

Since only the ventricles were visible in our images, we discarded the atria and big vessels from each segmented shape.  Each resulting mesh has 4446 vertices (the has 3052 vertices and the has 1776 vertices, the right-most part of the septum belongs to both ventricles) in point-to-point correspondence, and 9004 triangles.  Only the frame is selected for the analysis.  Since meshes are in point-to-point correspondence and share the same connectivity, they can be analysed using the [@Cootes1995a].  In , the shape of each patient $j$ is associated to a vector with the concatenated position $(x, y, z)$ of its nodes, giving a shape vector: $$X_j = ( x_0, y_0, z_0, x_1, y_1, z_1 \dots x_N, y_N, z_N)$$

We applied generalised Partial Procrustes Analysis [@Dryden1998a] to align the meshes thus removing the positioning and orientation variability.  We maintained size during this step, but differences in size due to anthropometric variables will be identified and removed in other steps of the framework.  This algorithm is an iterative method, at each iteration computes an estimation of the mean shape $\bar{X}^j$ and then rigidly registers each shape to that estimated mean.  This is repeated until convergence, in order to obtain an unbiased mean.  Here we show a full iteration $j$ of the algorithm: $$\bar{X}^j = \frac{1}{n} \sum^n X_i^j$$ $$\forall i: \; R_i^j, t_i^j = \operatorname*{arg\,min}_{R, t} \|R (X_i - t) - \bar{X}^j \|^2$$ $$X_i^{j + 1} = R_i^j(X_i - t_i)$$

Confounding deflation {#methodology:subsection:residual}
---------------------

To identify and remove the shape variability related to the confounding variables ($M$), and not to the studied condition, we use a procedure similar to partial correlation, as done in [@Singh2014a] for a regression problem.  We assume that shape $X$ can be decomposed in the sum of the population mean ($\mu_X$) and some deformation from that mean that is composed of noise ($\epsilon$), variability caused by the confounding variables ($X_M$) and variability from other sources ($X_i$).  The last component includes remodelling due to pathology : $$X = \mu_X + X_M + X_i +  \epsilon$$

The first step consists in estimating $X_M$ by building a linear model (whose coefficients are $w_M$) that predicts the expected shape from the confounding variables $M$. The training data of this model is only one population (the controls), to avoid introducing possible inter-population differences in the model.  This prediction will be our $X_M$. Then, for the shape of each individual, we subtract the predicted shape from the actual shape obtaining the prediction residuals.  Residuals represent the part of the shape variability that cannot be explained by the confounding variables.  To maintain the residual vectors in the same range as the original shape vectors, we summate the original population shape mean.  The final formula reads:

$$X_{res} =  X - X_M =  (X - w_c\cdot M) +\mu_X$$

The regression model coefficients $w_c$ are computed using with the shapes $X$ and confounding variables $M$. is a regression method that projects the input and output data to two low-dimension subspaces (called embeddings) that have maximal covariance [@Wegelin2000].  The confounding variables are standardised to have $1$ standard deviation and $0$ mean, but shapes are only standardised to have $0$ mean.  The full process is described below, and consists in an iterative process where at each iteration a new regression dimension of the low-dimensional spaces is computed.  The predicted part is then removed from the input and output spaces, and the process is repeated until the desired number of iterations is reached.  Several versions of exist, and we used Wold’s version.  All versions agree in the first iteration, but give different result for embedding spaces of more than one dimension.  Here we show the algorithm for the $r$-th iteration of a associating an input space $X$ with some response $Y$, both being matrices.  In our setting, and only for this part $X$ will be the confounding variables ($M$ in the rest of the paper) and $Y$ the shape vectors ($X$ in the rest of the paper). 

First, we compute the new dimensions of the embedding spaces by solving an eigenvector problem: $$\label{method:eq:dimRedPLS}
u^r, v^r  = \operatorname*{arg\,max}_{\|u^*\| = 1,  \|v^*\| = 1} {{u^*}^t X^t Y v^*}$$

We compute the rank one approximations (hat variables), and the regression coefficient $w_r$; for OLS refers to the classical ordinary least squares.  $$\begin{split}
\hat{X}^r = OLS_{prediction}(X^r u^r, Y^r v^r) \\
w^r = v^r (u^r)^t  \frac{(X^r u^r)^t Y^r v^r }{X^r u^r,} \\
 \hat{Y}^r = OLS_{prediction}(Y^r u^r, X^r v^r )
\end{split}$$

These rank one optimisations are used to update the input $X$ and its output $Y$ by removing the part that is already predicted.  $$\begin{split}
X^{r + 1} = X^r - \hat{X}^r  \\
Y^{r + 1} = Y^r - \hat{Y}^r \end{split}$$

The final result of the algorithm is the regression coefficients $w$ from $X$ to $Y$; and it is obtained by summing all the one-rank approximation regression coefficients.  $$w = \sum_r w^r$$

The coefficients $w_i$ of the prediction model, associated to confounder variable $M^i$, of this model can also be visualised and interpreted.  We can consider the partial determination coefficient, and visualise the shape pattern most associated to a certain confounder, as if it were a discriminating shape (see subsection \[methodology:subsection:classification\]). 

### Dimensionality reduction

Given the high dimensionality of the shape vectors and low number of samples, we use a method to find a subspace that contains the most relevant shape patterns.  The general linear model reads:

$$X_{orig} = \mu_X  + K X_{red} + \epsilon$$

Each shape $X_{orig}$ is expressed as the population mean($\mu_x$) plus some shape-specific associated low-dimensionality vector $X_{red}$ and a noise term $\epsilon$. The embedding matrix $K$ is constructed depending on the dimensionality reduction method, and contains the most interesting shape directions according to a certain metric. 

We use three different methods of linear : , and a combination of both.  and have been reported to be used in conjunction to classification methods, and in particular logistic regression [@Bastien2005].  We have described the regression modality of in subsection \[methodology:subsection:residual\], but also computes a space.  The difference is that we discard the regression coefficients $w$ and only use embedding vectors $u^r$ in equation \[method:eq:dimRedPLS\] of the input space.  These vectors define a vector subspace, but are not guaranteed to be orthogonal, so we use $QR$ decomposition [@Lord2007] to obtain an orthonormal base.  The combined method is based on prefiltering the shape using , keeping a high number of components ($>90\%$ of the variance), to then use as input to a .  This decreases the computation time, and adds stability by denoising the data.  Contrary to the typical procedure in machine learning, we chose not to standardise the modes by variance before applying , as the variance of each mode carries important information of the signal-to-noise-ratio. 

Classification {#methodology:subsection:classification}
--------------

In this section, we present the method used to train the classifier model and use the model coefficients to obtain the most discriminating shape.  The shape features obtained from the are combined with the confounding variables in a logistic regression model.  We choose logistic regression because we expect not to have a complete separation between both populations, and we want the model to be simple and interpretable.  The logistic model gives a probability that an individual $j$ with shape $X_j$ and confounding variables $M_j$ belongs to the pathological or control populations.  $$Pr(y_j = \textit{control}\ | X_j, M_j) = logit(\langle X_j, w_X \rangle + \langle M_j  w_M \rangle + b)$$

Logit refers to the logistic function $x \mapsto 1/(1 + exp(-x) )$.  $w_X$ and $w_M$ are the logistic regression coefficients for the shape and confounding variables respectively, and are chosen to minimize the log-loss of the probability of the training data $X$. The log-loss is the logarithm of the probability that the model is inconsistent with the observed data:

$$\label{eq:logloss}
\frac{1}{n}\sum^n \log \left(\left| y_j - Pr(y_j = 1 | X_j, M_j) \right|\right)$$

The logistic regression coefficients associated to the shape $w_X$, can be mapped from the reduced shape space to the full space by using the pseudoinverse of the dimensionality reduction matrix.  Let $w_{X_{red}}$ be the coefficients associated to the reduced shape models, and $K_{PCA}$ and $K_{PLS}$ the projection matrices of and respectively.  Then the coefficients associated with the full shape are: $$w_X = K_{PCA}^t K_{PLS}^t w_{X_{red}}$$

Then, we can visualise and interpret the shape pattern.  Analogously to multivariate regression, where the coefficients are indexed by the standard deviation to allow comparison among them, we need to adjust for differences in variance of the different coordinates.  Since node coordinates carry no meaning on their own, we treat shape as an object itself and do whitening of the shape, as is typically done with other multidimensional signals [@Kessy2018].  Since we are only interested in the remodelling direction, we normalise the vector to be unitary in the $L_2$ norm.  With these corrections, all the shape features are correctly scaled by their importance in prediction.  The full process to find the standardised shape pattern $\hat{w}$ reads:

$$\hat{w} =\frac{ \Sigma^{1/2} w_X}{ \| \Sigma^{1/2} w_X \|}$$

Where $\Sigma$ is the covariance matrix estimated using .  For visualisation, we can generate shapes that are representative of that shape pattern by adding the remodelling shape pattern, scaled with a parameter $\lambda$, to the mean shape.  To keep the shapes within the original range, we impose that $\lambda$ has to be within $3$ standard deviations of the variance associated to the shape pattern. 

$$\label{method:eq:syntheticMeshes}
X_{repr} (\lambda) = \mu_{X} +  \lambda \hat{w}$$

We can quantify the presence of remodelling in each shape, obtaining a scalar score for each individual, by computing the dot product of the shape vectors with the raw logistic regression coefficients associated to the shape only.  The previous whitening is only done for visualisation and comparison of modes.  If the shape pattern needs to be quantified in a population, the original one without standardization needs to be used. 

$$\textit{score}_i = \langle w_X, X_i - \mu_X  \rangle$$

The shape patterns can be compared using the $L_2$ dot product between standardised shape vectors, which coincides with the correlation of the scores associated to each pattern. 

Experimental setup
==================

Clinical dataset
----------------

The study comprises 77 controls and 89 athletes that underwent a , to study the cardiac remodelling triggered by the practice of endurance sport.  The study was approved by a local ethical board, and all participants gave written informed consent for the handling of their data.  Recruited athletes had been training an endurance sport, triathlon, over 10h a week during the last 5 years.  Table \[clinical:demographics\] shows the demographics of both populations. 

The controls and athletes come from different studies, and the demographics of both populations did not match exactly in age, but roughly represent the same general population in age and gender.  The study protocol and radiologist were the same for both cases.  Age is statistically different between athletes and controls, but the difference is very small (2 years).  We do not expect big differences due to this imbalance, since both athletes and controls are middle-aged.  There are also statistically significant differences in both weight and , but these correspond to physiological remodelling since endurance athletes are obviously fitter than the general population.  We used as possible confounders age, and gender, the typical adjustment variables in cardiology studies. 

The acquisition was ECG-gated from the R-peak during breath-hold.  The machines were Siemens Aereo and Siemens Magnetom, with an in-plane spatial resolution ranging from 0.5mm to 1mm.  The spacing between slices range between 8mm and 10.4mm, and the slice thickness was 8mm.  sequences were acquired with 25 frames per cardiac cycle.  The ventricular contours (epicardium and endocardium in the case of the , and only epicardium of the ) were automatically segmented from the SA using the automatic procedure described in the methodology. 

The quality of the automatic segmentations was assessed by one of our experts, but no manual refinement was performed in order to preserve point-wise correspondences.  Cases where errors could not be considered to be small were discarded: two individuals (both of them athletes) were discarded because the segmentation was inconsistent with the image.  The segmentation and registration errors are handled as noise in our study.  The meshes were very uniform and we found no self-intersecting artefacts.  As a consequence of the thick slices, the apex was not correctly imaged and presented much more noise than the basal part of the ventricles. 

                               Athletes     Controls     p-value
  ---------------------------- ------------ ------------ -----------
  Age [\[]{}y[\]]{}            35.4(6.1)    33.4(3.8)    0.013
  BSA [\[]{}m$^2$[\]]{}        1.78(0.19)   1.86(0.20)   0.005
  Weight [\[]{}kg[\]]{}        66.8(11.3)   73.5(15.1)   0.001
  Height [\[]{}m[\]]{}         1.71(0.09)   1.73(0.08)   0.151
  Women [\[]{}%[\]]{}          0.48         0.44         0.938
  Resting HR [\[]{}bpm[\]]{}   57.2(8.4)    65.8(10.6)   &lt;0.001

Automatic measurements
----------------------

We computed automatic measurements of the 3D shapes that are analogous to the classical clinical measurements, using the point-to-point correspondence and labelling coming from the model-based registration and segmentation.  This allowed a better understanding of the discriminative shape patterns by assessing how these measurements vary in response to the remodelling score $\lambda$ on the synthetically generated meshes according to equation \[method:eq:syntheticMeshes\].  We computed the volumes of both ventricles, as well as the myocardial mass of the . 

BMI-based downsampling
----------------------

Obesity (defined as &gt; 30) and overweight (defined as &gt; 25) have been reported as risk factors to cardiovascular illnesses in the literature and have a clear influence in cardiac shape and function [@Alpert2018].  Surprisingly, overweight and athletic remodelling share similarities.  Even if one might expect them to be opposite, as endurance athletes and overweight body fat are in the opposite sides of the spectrum, both remodellings are triggered by an increase of the heart’s loading.  In the case of the athletes, it is the increase of during exercise that produces a volume overload, and in the obese it is a mix of increased needs at rest to account for the bigger body size [@Lavie2007] and a ventricular pressure overload due to an increase of arterial pressure [@Messerli1982]. 

To study the effect of an imbalance in a shape-affecting variable between the control and case populations, we biased our population to increase the of the controls, as overweight has a well known effect on the heart.  To generate the imbalance, the control class was downsampled to maintain only 25% of its controls, favouring keeping the ones with higher .  The individuals to remove were selected randomly among the controls, with a probability of being kept proportional to the rank of its in the control population.  Athletes were not downsampled.  This procedure was repeated with 100 different seeds to obtain different imbalanced datasets and add robustness to the results. 

We analysed how this imbalance affected the remodelling pattern found, and to which extent could confounding adjustment and confounding deflation correct this effect.  To study the stability of the shape pattern, we computed the $L_2$ product between the discriminative shape patterns obtained with the downsampled datasets and the results obtained with the full dataset, that serves as groundtruth.  This was performed for the different methods.  We also tested the $L_2$ product of the discriminative patterns with the -shape pattern, obtained with an adaptation of our framework to the regression problem.  Additionally, we compared qualitatively the automatic measurement response to the remodelling score of both the downsampled and original dataset.  We tested both the covariate adjustment, and the confounding deflation.  For the confounding deflation, we also evaluated how different choices during the training of the shape prediction model affected the obtained most discriminative shape pattern. 

Results
=======

Dimensionality reduction
------------------------

To choose the best configuration of parameters for the method, we used 10-fold and computed the mean log-loss of all the validation set, defined in equation \[eq:logloss\], over a wide combination of parameters.  In table \[experiments:table:cvTableReduced\] we find the log-loss of the best parameter choice for the 3 different methods (PCA, PLS and PCA + PLS).  They correspond to a PCA with 5 modes, a PLS with 3 and the PCA + PLS with 20 modes and 3 modes.  This experiment also provided an overview of how the use of demographics affected the classification metrics: there is a considerable improvement in terms of log-loss when the confounding variables are used in the model, and a minor improvement when using PLS instead of PCA.  When using PCA + PLS, the metric is very similar to PLS.  Confounding deflation gave a worse result than the raw (non-deflated) shapes when adjustment is used.  However, confounding deflation improved the resulting classification metric with respect to the raw when confounders are not added to the logistic regression model. 

  ---------------------- --------- ------------------ --------- ------------------
                                                                
  Method                   No adj.    Confounders adj.    No adj.    Confounders adj. 
  PCA$_5$                     0.58               0.46      0.48               0.46
  PLS$_3$                     0.59               0.44      0.52               0.48
  PCA$_{20}$ + PLS$_3$        0.59               0.43      0.50               0.46
  ---------------------- --------- ------------------ --------- ------------------

T

Athletic model 
---------------

We applied our framework to identify the athletic remodelling in our dataset, and compared the effect of using the different method and confounding-bias correction methods (confounding deflation and adjustment).  We used the parameters found in the previous section via .  Figure \[experiments:figure:correlationAllModels\] shows the $L_2$ product between all combinations of methods and the adjustment or not by confounding variables.  For the original (no confounding deflated) shapes (left), we found big differences due the inclusion or not of the confounding variables.  The choice gave only minor differences in this dataset: the resulting most discriminating modes of the were very similar.  When the confounding deflation was applied (Figure \[experiments:figure:correlationResidualModels\]), the use of the confounding adjustment became redundant and did not influence the resulting shape. 

[0.4]{}

[0.4]{}

Figure \[methodology:figure:fullPopulationPCA\] depicts the most discriminative shape patterns, expressed as mean + 2 STD, obtained by the model with and without adjustment.  Figure \[methodology:figure:fullPopulationMeasurements\] shows their asociated measurement response.  Both the adjusted and unadjusted models found an increase of the ventricular volumes that is the same for both and .  Both models found shape changes in the : the outflow region dilated more than the inlet base and apex.  Regarding myocardial mass, however, the models gave different results.  After confounder adjustment, we observed a large increase of myocardial mass similar to ventricular dilation, while the unadjusted model found no difference in the mass.  The increase of mass of the adjusted model is concentrated in the base and lower in the apex: the base of the is flatter and is exposed to bigger mechanical stress, so it is more prone to compensatory myocardial thickening.  Finally, the increase of ventricular volume in the adjusted model is more pronounced than in the unadjusted model. 

The adjustment or not by confounders resulted in big differences between the remodelling patterns found.  On the other hand, differences due to were smaller, yet noteworthy.  Figure \[methodology:figure:differenceDRModels\] shows the remodelling patterns found for the different techniques, with a colormap showing the local amount of remodelling: reed indicates substantial changes and blue no remodelling.  The figure shows a view of the free wall, which is the region that experienced more shape changes.  All 3 shape patterns were similar and followed the same global trends, but they presented regional differences.  We observed that both and + found a remodelling that was localised in the outlet and , at a lower degree, in the apex, while showed a more spatially distributed remodelling that also affected the base and had a smoother transition between the affected and unaffected areas. 

BMI model {#experiments:regressionBMI}
---------

To have a better understanding of the effect of an elevated in the ventricles, we constructed a regression model that predicts from the cardiac shape.  As athletes and controls had different ranges of , and to avoid finding any interference with the remodelling due to endurance sport, we only use the controls to build the model.  Also, since and body size are related, we did not use the confounding variables in this model.  The model was built using a in its original regression mode, using $3$ dimensions.  To evaluate the model prediction capability, we computed the determination coefficient $R^2$ using 5-fold , which was $0.44$. Similar to classification, we obtained synthetic representative shapes ($\hat{x}_b$) associated with certain value $b$. The synthetic mesh associated with a certain was the one having minimal-distance to the mean shape (using the Mahalanobis distance), constrained to having the required predicted . 

The representative shapes associated to BMI values of 17.5 and 30, which are very extreme values, can be seen in Figure \[experiments:figure:bmiModel\], and Figure \[experiments:figure:bmiMeasurements\] shows the measurement response of the remodelling.  The -associated remodelling consists of a moderate increase of ventricular size and a bigger increment of the myocardial mass. 

Confounding adjustment
----------------------

Figure \[experiments:figure:stabilityRaw\] depicts the $L_2$ products between the athletic remodelling shape pattern (normalised to be a unitary vector) derived from the full population, and the shape pattern derived from the -imbalanced populations.  We found that without adjustment the imbalance confused the method, and it mixed the differences due to with the ones related to endurance training.  The adjustment by confounding variables resulted in a better agreement between the downsampled-derived remodelling patterns and the full-population pattern than the unadjusted model, as can be seen in Figure \[experiments:figure:measurementsConfRaw\].  Figure \[experiments:figure:bmiRaw\] shows the dot product between the most discriminating shape pattern obtained with the downsampled (boxplot) and full (thick solid lines) datasets, and also the shape pattern associated to , computed as described in section \[experiments:regressionBMI\].  There, we can see that in the downsampled datasets, the -associated remodelling and the athletic one had lower correlation.  This drop of correlation was larger for the models unadjusted by confounder variables. 

Figure \[experiments:figure:measurementsConfRaw\] shows the measurement response to the athletic remodelling of several models trained with a downsampled population.  It shows different methods with confounding adjustment (upper row) and without confounding adjustment (lower row).  We can see that when no confounding correction is used, the found is associated with lower myocardial mass, contrary to what is observed with the full dataset experiment in which myocardial mass was maintained.  When the shape is adjusted by confounders, athletic remodelling was again associated with an increase of myocardial mass, consistent with the full population model.  Figure \[experiments:figure:athletesModelsDownsampled\] shows the discriminating shape patterns for a randomly selected downsampled dataset.  There, we can observe that the model associates athletic remodelling with a decrease of the wall thickness in the septum and apical regions of the when confounders are not considered. 

Resulting shape patterns were affected by the methods, was more unstable than and + , who had a better correlation to their full-population discriminating shape.  In Figure \[experiments:figure:measurementsConfRaw\], we can see that was not able to recover the increase of myocardial mass, and its mode presented no variation in the myocardial mass while both and + could. 

[0.45]{}

[0.45]{}

[0.8]{}

\

[0.8]{}

Confounding deflation {#confounding-deflation}
---------------------

In this subsection we analyze the effect of the confounding deflation.  As stated above, the shape prediction model of the confounding deflation was trained using only the controls.  We tested two possible scenarios: the case in which the athletes population was downsampled, and therefore the building of the shape residual model was not affected; and when the controls were downsampled (as in the previous subsection).  The former was the most appropriate situation to apply confounding deflation, since there would not be extra unstability/bias introduced during the confounding deflation step, while the latter could introduce the bias in the dataset during the confounding deflation process.  Finally, we show the potential danger of training the residual model using both populations and how the confounding deflation could even increment bias. 

When controls were not downsampled, the population used to train the shape prediction model was relatively unbiased.  Instead of downsampling the controls, we downsampled the athletes, removing athletes with high analogously to the procedure used to downsample the controls.  Figure \[experiments:figure:stabilityRes\] and \[experiments:figure:bmiRes\] depict the same experiments as in the previous section: the dot product between the most discriminating shape obtained with the downsampled data and the remodelling obtained considering the whole population, and also the dot product with the mode.  Results showed a considerable decrease in variability and a higher correlation with the full-dataset-derived remodelling than the confounder adjustment experiments (Figure \[experiments:figure:stabilityRaw\]).  Adding confounder adjustment on top of confounding deflation did not produce an increase of accuracy. 

Figure \[experiments:figure:stabilityResDownControls\] shows the results when the population in which the residual is trained is downsampled.  We downsampled the controls based on their .  Results were much worse than when athletes were downsampled, and even worse than a simple confounder adjustment.  Adding confounding adjustment on top of confounding deflation had a beneficial effect.  Figure \[experiments:figure:stabilityResDownsampledBothPopulations\] shows the effect of using both athletes and controls in the training of the shape-prediction model for confounding deflation: there was a drop in stability compared to the use of a single population when and + were used (Figure \[experiments:figure:stabilityRaw\]).  Strangely there was an improvement compared to the baseline (where confounders are completely ignored). 

Therefore we observed that using both populations in the training of the shape prediction model during the confounding deflation step resulted in worse results.  We explored the reason why using both populations can create a confounding effect, associating an imbalance in a variable to the inter-class shape differences.  This can happen even when the variable is not associated to any shape remodelling.  To illustrate this effect, we created a dummy synthetic variable which is just the athlete label plus Gaussian noise.  Obviously, by generation we knew that this variable did not have any direct relationship to shape.  To evaluate this effect, we computed the $L_2$ dot product with the most discriminating shape between athletes and controls and the shape pattern associated to the dummy variable in the shape-prediction model.  We repeated this full variable generation and $L_2$ products computation process $100$ times to remove randomness of the analysis.  Figure \[experiments:figure:artificialData\] (left) shows the distribution of this dummy variable for a certain seed.  We constructed two residual models one only with the controls, and the other with both populations.  Figure \[experiments:figure:artificialData\] (right) shows the distribution of the dot product of the shape prediction model coefficients associated to the dummy variable when both populations are used, and when only the controls are used.  The shape associated to the dummy variable is independent to the athletic remodelling when only one population is used, but becomes very similar to the control-athlete difference when both populations are used in training due to confounding effect. 

[0.4]{}

[0.4]{}

[0.4]{}

[0.4]{}

Discussion
==========

The shape models corresponding to the athletic and overweight remodelling corroborated the current clinical literature.  The remodelling consisted in mostly a concentric remodelling, both the mass and EDV increased, increasing more the mass than the volume: this coincides with the hypothesis that remodelling is concentric to cope with elevated pressure (increase of myocardial mass) and demands (increase of volume).  Another hint that the remodelling was mostly pressure driven is that the septum of the mesh representing a high individual presents a bulge below the aorta, that has been described as an early indicator of elevated pressure [@Baltabaeva2008; @Gaudron2016].  The athletic remodelling consisted of a predominant increase of volume, and if confounder adjustment was performed, an increase of myocardial mass.  There is some controversy if endurance athletes remodelling is more eccentric ( they increase their volume more than their mass) or concentric (they increase their mass more than their volume), but there is a consensus that there is both an increase of mass and volume [@Scharhag2002], and therefore we considered the adjusted model more accurate.  The increase of ventricular volume in the confounder adjusted model is considerably larger than the unadjusted.  These differences between the adjusted and not adjusted models can be explained because the variability in size and myocardial mass can come also from the patients morphometrics: big persons have big hearts.  Without considering this extra non-imaging information, it is not possible to discern when the size of the heart is due to the patient being big, or the heart had a dilation as a remodelling reaction to exercise. 

The artificially generated imbalanced sets allowed us to validate our hypothesis: without any correction, the athletic remodelling presented a reduction in myocardial mass.  This is obviously false, since all studies have found that endurance exercise provokes an increase of myocardial mass.  The downsampled control population had a high percentage of overweight, who had a concentric remodelling and causing the previous association with controls and higher myocardial mass in the unadjusted model.  The stability analysis, also showed a bigger match between the downsampled-derived remodelling and the full-dataset one when confounding adjustment was used. 

We also studied another strategy: confounding deflation.  Confounding deflation consists of generating a model that predicts shape from the confounders, and working with the residual of that prediction.  Our results showed that this strategy worked well when there was access to a good population to train the shape prediction model: otherwise it can actually increase bias.  This was shown when training on the downsampled population, or when both the case and control were used simultaneously. 

Finally, we compared different linear methods for : and .  In the original dataset, we obtained a better classification accuracy using .  We also showed that was able to capture more localised remodelling than , that was limited to smooth and global remodelling.  However, in stability tested during the downsampling analysis.  Both models could be combined by using a coarse , that only removes the modes encoding very little variance, with a .  With this method we were able to capture both localised remodelling and have high stability to imbalances. 

Conclusion
==========

We have presented a framework to find regional shape differences between two populations, taking special care to correct for any potential bias related to demographics parameters.  The framework is fully linear: it used a sequentially combined with a as dimensionality reduction of the shapes, followed by a logistic regression model.  The linearity allows to easily interpret and visualise the model and build synthetic representative shapes of the model.  To correct for confounding effects, it incorporates adjustment, where the confounding variables are added to the logistic model and confounding deflation, which consists of building a regression model that predicts shape from confounding variables, and we used it to remove the shape variability associated to these variables. 

We applied our framework to a real dataset consisting of athletes and controls to find the remodelling due to the practice of endurance exercise.  Our results confirmed the current literature on endurance-sport remodelling in the : ventricular dilation and increment of myocardial mass, specially in the basal area.  In the , we found that the volume increase was not homogeneous but concentrated in the outflow.  In the controls, we used an adaptation of our classification framework to regression to explore obesity remodelling and found it to be mainly an increase of myocardial mass. 

In this population, we analysed the effect of confounders in a semi-synthetic dataset obtained by downsampling the control population non-uniformly, keeping individuals with high .  Even if athletic remodelling is very prominent, we were able to bias the model to output that athletes have lower myocardial mass than controls.  This was corrected when adjustment was used.  However, we found that we could only use confounding deflation when the control population was relatively big and balanced, and if that is not the case using confounding deflation can actually increment bias. 

In our work we have tested only linear methods, but this confounder-related problems might appear even more with the use of more complex frameworks, able to capture non-linear shape patterns to capture subtler morphology differences in populations that are almost indistinguishable from controls and differences cannot be found through traditional means. 

Acknowledgements {#acknowledgements .unnumbered}
================

We thank Dr Weese and Dr Groth from Philips Research Hamburg for the segmentation tool

Funding {#funding .unnumbered}
=======

This study was partially supported by the Spanish Ministry of Economy and Competitiveness (grant DEP2013-44923-P, TIN2014-52923-R; Maria de Maeztu Units of Excellence Programme - MDM-2015-0502), el Fondo Europeo de Desarrollo Regional (FEDER) , the European Union under the Horizon 2020 Programme for Research, Innovation (grant agreement No.  642676 CardioFunXion) and Erasmus+ Programme (Framework Agreement number: 2013-0040), “la Caixa” Foundation (LCF/PR/GN14/10270005, LCF/PR/GN18/10310003), Instituto de Salud Carlos III (PI14/00226, PI17/00675) integrated in the “Plan Nacional I+D+I” and AGAUR 2017 SGR grant nº 1531.