In medical imaging almost every image is at least weakly labeled with clinical indicators. A minimal set of such labels usually includes age, sex, demographics, and diagnostic or prognostic clinical values such as blood pressure, history of smoking or outcome scores. The labels are of interest for clinical decisions or clinical research, however, they are rarely included in the analysis of shape or image features. We propose to construct joint statistical shape and attribute models to better understand and visualize the interplay between shape deformations and the variations of clinical indicators. The proposed model promises to reveal associations between anatomical shapes, image features, and clinical phenotypes, to help understand disease mechanisms and risk factors, and to improve predictions and patient stratification in clinical trials.
We demonstrate our approach in the context of a clinical cohort of ischemic stroke patients, where we investigate the relationships between the patterns of expansion of the lateral ventricles, the spatial patterns of white matter hyperintensity (WMH) indicative of the white matter disease, and patient data that include stroke outcome scores. Our model is built from segmentations of ventricles and white matter hyperintensities in T2-FLAIR magnetic resonance images.
We integrate copula  and Gaussian processes  to jointly analyze shape and other variables. Copula methods provide a decomposition of every continuous and multivariate distribution into its marginal distributions and a copula that captures the dependency structure 
. The separate modeling of the marginal distributions and of the copula offers better flexibility than the classical statistical shape models. The marginal distributions can be modeled as Gaussian, or other parametric and non-parametric distributions. The choice of the model for each marginal distribution is separate from all others, e.g., a Gaussian distribution for the shape components and the empirical distributions for the clinical indicators. In this work, we use a Gaussian copula, which means we transform the data into a space where all marginal distributions are Gaussian. In this latent space, we then apply standard shape modeling techniques, such as the principal component analysis (PCA) or the Gaussian processes, to model the dependency structure between 3D coordinates of the shape and other indicators. The transformation to the latent space is invertible and our shape model remains generative. We can, therefore, construct the conditional distribution given clinical indicators using Gaussian process regression. Here we focus on the setting where some clinical indicators are given and we capture the conditional distribution of shapes, i.e., a patient-specific statistical shape model.
The main contributions of this paper are: (i) a novel approach for generating conditional shape models based on clinical indicators or image features; (ii) copula-based shape models that include binary, discrete, ordinal and continuous variables; (iii) novel associations between ventricle growth, WMH patterns and clinical indicators in the ischemic stroke patient cohort.
Related Work Although many clinical indicators are collected and available with medical images, the idea of combining statistical shape models with clinical indicators is rarely explored. The work of Blanc et al.  includes patient attributes in femur modeling. The method assumes the data is Gaussian distributed and cannot handle binary, discrete, or ordinal variables. Pereanez et al.  proposed constraining shape models based on clinical indicators. They demonstrated that such constrained models are beneficial for model-based image segmentation. This approach is limited to a model prediction given a complete set of clinical indicators and handles ordinal variables as a set of binary attributes. Our method combines these ideas to employ the copula and uses the conditional distribution to capture the remaining variability given partial data. Han and Liu 
proposed a variant of PCA called Copula Component Analysis that is inherently scale-invariant, more robust to outliers and handles arbitrary discrete ordinal marginal distributions. Hoff
presented the idea of an extended rank likelihood which can also account for non-continuous variables. In computer vision, Eggeret al. [4, 5] used copulas for shape and appearance modeling to account for the non-Gaussian appearance of human faces, however, they did not handle discrete or ordinal variables.
Conditional shape models without clinical indicators are well studied. A Bayesian approach was proposed by Albrecht et al.  to model the remaining variability of the femur given partial information. They showed that (probabilistic) PCA and Gaussian process regression with a statistical kernel are equivalent for such modeling. We build upon this model and include clinical indicators to construct patient-specific conditional shape models based on image features and/or clinical indicators. In our experiments, we visualize the conditional models and demonstrate superior performance in a reconstruction task when including all data in the statistical model.
We propose a joint model for anatomical shape, image features and clinical indicators based on the copula  and combine it with the idea of conditional shape models based on Gaussian processes . To enable the use of not only continuous, but also binary, discrete and ordinal variables, we use a sampling strategy similar to .
We define an instance of our model to be a combination of 3D shape coordinates and image-based feature associated with surface point . All surface points are in correspondence across training data, achieved via shape registration to an atlas in our experiments. In addition each instance includes binary, discrete, ordinal or continuous variables . Let
be the vector representation of the instance,with . The training set of shapes with image features and clinical indicators forms the data matrix
. Let random variablerepresent component of the instance and represents its marginal distribution.
provides the dependency structure in the joint cumulative distribution function (CDF)such that
Where is the inverse standard normal CDF and is the joint CDF of a zero-mean multivariate Gaussian distribution with covariance matrix . The multivariate latent representation
distributed according to the standard normal distribution.
To learn the model above from training data we first estimate the marginal distributions with parametric or nonparametric methods. For properties that cannot be naturally captured by a Gaussian distribution (e.g., image features, sex or outcome scores), we employ the empirical marginal distribution for all the components as proposed in
. If only limited data is available, various assumptions on the marginal distribution must be made, separately for each component. In statistical shape modeling, a Gaussian marginal distribution is usually assumed. To transform our data into the latent space, we first transform the empirical marginal distribution into a uniform distribution by replacing every entry of the data by the CDF of the empirical distribution. Second, we use the inverse CDF of a standard normal distribution to map to the latent space.
While the original copula framework is limited to continuous distributions, multiple solutions have been proposed to overcome this limitation. The first step, the transformation to a uniform distribution, relies on unique ranks (values of CDF). For binary, discrete and ordinal variables this sorting is not unique. In our experiments, we observed that the latent representation is stable when permuting elements with the same values for ranking, in the setting where we have many more variables than samples (). To construct a consistent dependency structure we generate multiple (50) random rankings for non-unique values and average the resulting estimates of the model parameters. Once the data is transformed to its latent representation we estimate the covariance matrix and represent it by its eigenvectors (i.e. principal components).
Conditional modeling The Gaussian copula represents the dependency structure of our data based on the latent representation . We use Gaussian processes to capture correlations in the latent space . We define function on a finite domain such that refers to the th element of the th training example in our latent space, , . The Gaussian process is fully defined by the mean and covariance functions estimated from the training data:
Given a subset of the vector as an observation, we construct component-wise its latent representation . The prediction of the missing parts reduces to a standard Gaussian process regression problem. We assume observations are corrupted by uncorrelated Gaussian noise and obtain the prediction for the distribution of the complete vector from the observed vector . The conditional distribution is again a Gaussian process with mean:
is the identity matrix.
The parameter has an impact on the accuracy of the prediction and varies between different observations. The different modalities in our model have a different amount of uncertainty. For shape observations, is an estimate of the segmentation accuracy and the registration quality, for clinical indicators, it is a measurement of how accurate they are measured. We estimated the optimal parameters for observation uncertainty for each set of observations via cross-validation.
We produce a statistical model based on a dataset of T2-FLAIR magnetic resonance images from 793 acute ischemic stroke patients, without gross pathologies, such as mass effects, from the Genes Affecting Stroke Risk and Outcomes Study (GASROS), for which patients were enrolled at Massachusetts General Hospital between 2003 and 2011 . Manual WMH segmentation was performed using MRIcro , and ventricles were automatically segmented (Dice: 0.89), using a 3D U-Net . Subsequently, we isolated periventricular from subcortical WMH, based on connectivity of individual WMH lesions to the ventricles. After segmentation, we establish correspondence between every ventricle segmentation and the atlas ventricles by fitting a Gaussian process deformation model estimated via the iterative closest point algorithm . The ventricles are represented as a surface triangle mesh with 1504 vertices. The WMH is modeled as a feature based on the vertices of the mesh. Every voxel labeled as periventricular WMH is assigned to the closest vertex point. This results in a WMH representation as a surface feature. In addition to the shape of the ventricles and the WMH burden feature, we include clinical indicators that represent age, sex, hypertension (binary), hyperlipidemia (binary), atrial fibrillation (binary), smoking (ordinal, 1-6) and NIHSS (stroke score, ordinal 0-42), mRS (outcome score, ordinal 0-6) and the overall WMH volume. With 1504 3D ventricle vertices, 1504 WMH burden features and 9 clinical indicators, the total dimensionality of a single data point is .
Fig. 1 presents the joint model of the ventricle shape, WMH burden features, and clinical indicators and also three conditional models for different clinical indicators, like age and sex. The conditional models have similar mean and variability as the unconditioned model when conditioned on features with low correlations to the remaining components Fig. 1b and show different means and less variability when conditioned on stronger correlated Fig. 1c or combined features Fig. 1d. The supplementary material contains a video111https://youtu.be/gPoHP_iFQIA of our graphical user interface for data exploration. It enables the user to fix clinical indicators and visualize the remaining variability in the conditional model.
|mean||ventricles||WMH||indicators||ind + vol||combined|
|Age||15.77 9.49||11.27 6.62||12.06 7.58||13.21 8.19||11.57 7.04||9.85 5.76|
|mRS||2.01 1.50||1.88 1.37||1.96 1.45||1.65 1.20||1.65 1.20||1.65 1.19|
Clinical indicator prediction. Mean error and standard deviation are reported for the prediction of age, stroke outcome (mRS) and sex. The columns of the matrix correspond to the prediction based on the mean of the marginal distribution and conditioned on the ventricle shape, WMH burden, all other clinical indicators, indicators combined with WMH volume (ind + vol), and all available components except for the one being predicted. For sex, the percent of correct predictions is reported.
|mean||ventricles||WMH||indicators||ind + vol||combined|
|ventricles||1582.2 34.5||-||132.9 35.0||146.9 45.1||144.8 43.0||132.3 34.8|
|WMH||3980.1 2830||2927 2270||-||3073 2419||1291 897||1254 866|
We evaluate the joint model in a reconstruction task of estimating the full vector from partial data in three different tasks: attribute prediction (Table 1), shape prediction, and WMH burden prediction (Table 2). For all tasks, we use the mean of the conditional distribution in the latent space to generate predictions. The values are obtained by randomly splitting the dataset into a training set of 600 patients and a validation set of the remaining 193 patients. The attribute prediction task (Table 1) aims to estimate clinical indicators, e.g., stroke outcome, given the shape and WMH burden. We observe that the ventricle shape and WMH burden are strong indicators for age whilst they contribute little to the stroke outcome score. For age and sex, using the combination of all modalities provides the best prediction. The shape and WMH prediction tasks (Table 2) reveal how the imaging data from a patient relates to the clinical indicators. The shape of the ventricles can be predicted from the WMH burden. The prediction in the other direction is more challenging, due to strong variability in WMH burden in the population. Clinical indicators including WMH volume combined with the ventricle shape produce the best prediction.
In the third experiment, we demonstrate the difference between the unconditional and conditional models. We produce the sample distributions by drawing 1,000 random samples from the corresponding model. Fig. 2 reports those sample distributions before and after conditioning on the same three events as in Fig. 1. We observe that the conditioning on male has a weak effect whilst the conditioning on age or a combination of age and sex has a strong effect on the resulting marginal distributions.
We implemented the copula model based on the code provided in  in the Scalismo222Scalismo - A Scalable Image Analysis and Shape Modeling Software Framework
https://github.com/unibas-gravis/scalismo framework for statistical shape modeling based on Gaussian processes. The computation of the conditional shape models based on the clinical indicators takes a few seconds for our dataset on a single core. Reducing the model via PCA will speed it up to real-time performance, enabling interactions with the model.
Limitations The proposed model has two main limitations. First, we assume dense correspondence on the surface of the ventricles as well as correspondence across image features. This is especially limiting if some shape features have weak or no correspondence across the population. Second, the Gaussian model for the dependency structure we are using is limited to second-order dependency and does not reflect higher-order dependency. These two limitations and assumptions are very common in the statistical shape modeling, especially when working with small datasets.
We present a joint generative model of shape, image features, and clinical indicators. Conditional models for partially observed data can be straightforwardly derived. The copula can be implemented as a simple pre- and postprocessing step in existing pipelines relying on (probabilistic) PCA or Gaussian processes. In addition, conditional models can be constructed efficiently for interactive data exploration. Our work provides a highly effective approach to exploring high dimensional data, especially spatial correlations, for example between ventricle expansion and WMH burden patterns. Statistical shape models are commonly used as statistical priors for image analysis tasks, and the proposed conditional models can be interpreted as patient-specific statistical shape models for patient-specific medical image analysis. Our model includes binary, discrete, ordinal and continuous variables and hence can handle a vast variety of available clinical indicators.
Aknowledgments This research was funded by SNSF P2BSP2_178643, NIH NIBIB NAC P41EB015902, NIH NINDS R01NS086905, Horizon2020 753896, De Drie Lichten 24/18, ZonMw 104003005, Wistron Corporation, AWS, SIP and NVIDIA.
-  Albrecht, T., Lüthi, M., Gerig, T., Vetter, T.: Posterior shape models. Medical image analysis 17(8), 959–973 (2013)
-  Blanc, R., Seiler, C., Székely, G., Nolte, L., Reyes, M.: Statistical model based shape prediction from a combination of direct observations and various surrogates: application to orthopaedic research. Medical image analysis 16(6), 1156–1166 (2012)
-  Dubost, F., de Bruijne, M., Nardin, M., Dalca, A.V., Donahue, K.L., et al.: Automated Image Registration Quality Assessment Utilizing Deep-learning based Ventricle Extraction in Clinical Data. arXiv e-prints arXiv:1907.00695 (Jul 2019)
-  Egger, B., Kaufmann, D., Schönborn, S., Roth, V., Vetter, T.: Copula Eigenfaces with Attributes: Semiparametric Principal Component Analysis for a Combined Color, Shape and Attribute Model, pp. 95–112. Springer (2017)
-  Egger, B., Kaufmann, D., Schönborn, S., Roth, V., Vetter, T.: Copula eigenfaces - semiparametric principal component analysis for facial appearance modeling. In: Proceedings of the 11th Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications: Volume 1: GRAPP. pp. 50–58. SCITEPRESS-Science and Technology Publications, Lda (2016)
-  Genest, C., Ghoudi, K., Rivest, L.P.: A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika 82(3), 543–552 (1995)
-  Han, F., Liu, H.: Semiparametric principal component analysis. In: Advances in Neural Information Processing Systems. pp. 171–179 (2012)
-  Hoff, P.D.: Extending the rank likelihood for semiparametric copula estimation. The Annals of Applied Statistics pp. 265–283 (2007)
-  Lüthi, M., Gerig, T., Jud, C., Vetter, T.: Gaussian process morphable models. IEEE PAMI 40(8), 1860–1873 (2018)
-  Pereañez, M., Lekadir, K., Albà, X., Medrano-Gracia, P., Young, A.A., Frangi, A.: Patient metadata-constrained shape models for cardiac image segmentation. In: STATCOM, MICCAI. pp. 98–107. Springer (2015)
Rasmussen, C.E.: Gaussian processes in machine learning. In: Summer School on Machine Learning. pp. 63–71. Springer (2003)
-  Rost, N.S., Rahman, R.M., Biffi, A., et al.: White matter hyperintensity volume is increased in small vessel stroke subtypes. Neurology 75(19), 1670–1677 (2010)
-  Sklar, M.: Fonctions de répartition à n dimensions et leurs marges. Université Paris 8 (1959)
-  Tsukahara, H.: Semiparametric estimation in copula models. Canadian Journal of Statistics 33(3), 357–375 (2005)
-  Zhang, C.R., Cloonan, L., Fitzpatrick, K.M., Kanakis, A.S., Ayres, A.M., Furie, K.L., Rosand, J., Rost, N.S.: Determinants of white matter hyperintensity burden differ at the extremes of ages of ischemic stroke onset. Journal of Stroke and Cerebrovascular Diseases 24(3), 649–654 (2015)