Randomized denoising autoencoders for smaller and efficient imaging based AD clinical trials

Our contributions

We achieve these goals by proposing a novel learning model based on deep learning. Deep architectures are non-parametric learning models [1,3] that have received much interest in machine learning and computer vision recently. Although powerful, it is well known that they require very large amounts of unsupervised data, which is infeasible in neuroimaging, where the dimensionality d of the data is always much larger than the number of instances (n). A naïve use of off-the-shelf deep learning models on neuroimaging data expectedly yields poor performance. In the last few months, however, independent of our work, deep learning methods have been used successfully in structural and functional neuroimaging [12,5,10]. To get around the difficulty highlighted above, [12] uses a region of interest approach whereas [5,10] sub-samples each data instance to increase n. Our work provides a mechanism where no such adjustments are necessary. The key contributions of this paper include (a) Scalable deep architecture(s) for learning problems in neuroimaging where number of data instances is much smaller than the data dimensionality (i.e., our models permit whole–brain analysis) and (b) An imaging derived continuous measure with smaller variance that leads to efficient AD clinical trials with moderate sample sizes (and based only on one time–point data).

2.1 Stacked Denoising Autoencoders (SDA)

We motivate our formulation by highlighting the difficulty in using stacked denoising autoencoders (SDA) [1,3] directly in the n ≪ d regime. An autoencoder is a single layer neural network that learns robust distributed representations of the input data. A denoising autoencoder (DA) constructs these representations by stochastically corrupting the inputs. Denoting the d dimensional inputs by ${\{{\mathbf{x}}_i\}}_1^n$, a DA outputs h_i = σ(Wx_i+p) (σ is a point–wise sigmoid) by minimizing the loss ℒ(․) between the input and its reconstruction x̂_i = σ(W^Th_i + q) as,

where γ(․) is the point–wise stochastic corruption [1]. A stacked denoising autoencoder (SDA) greedily concatenates L(> 1) DAs, i.e., l^th layer outputs are the un–corrupted inputs for (l + 1)^th layer,

where θ denotes the full set of stochastic gradient (SG) learning parameters (corruption rate, learning rate, hidden layer length). The transformations ${\{{\mathbf{W}}^l,p^l,q^l\}}_1^L$ serve as a warm–start for supervised tuning where one compares the output of the L^th layer to ${\{{\mathbf{y}}_i\}}_1^n$. This greedy layer–wise unsupervised training followed by supervised fitting is central to most deep architectures [3,1].

Recall that SG learning is expected to converge to a local minimum only in the asymptotic setting (of large n). Hence, the warm–start described above is only reliable when large amounts of unsupervised data are available, which is the case in computer vision but not in neuroimaging. Otherwise, the network overfits whenever d is much larger than n. In neuroimaging, d is generally on the order of millions (number of voxels) and n < 1000. Hence, traditional SDAs cannot be directly used (they will generalize poorly). Recent work uses deep architectures in neuroimaging either by reducing d (using anatomical ROIs or feature selective) or increased n (splitting a data instance using sets of 2D slices) [12,5,10]. Nonetheless, frameworks to perform whole–brain analysis (the de-facto input when SVMs are used for brain image classification) will yield improvements by exploiting 3D local neighborhood dependencies directly.

2.2 Randomized Denoising Autoencoders (rDA)

In Section 1, we motivated the task of concurrently optimizing two goals. Our system should be able to capture differences across different dementia stages (i.e., controls, MCI, AD) while at the same time keeping intra-stage prediction variance as small as possible. In other words, we seek to decrease the prediction variance at no cost of approximation bias. Although, these seem like competing requirements, it turns out that this behavior is exactly what is offered by ensemble learning [2]. Recall that Ensembles are bootstrap randomizations around sets of weak learners which reduce the prediction variance in expectation. So, properly incorporating an ensemble approach within a deep architecture should yield the behavior we expect. We can generate the ensembles for a given learner in multiple ways [2] — a randomization over the number of features and/or data instances. Here, we already have n ≪ d, so randomization over n is infeasible. Instead, we distribute/randomize over the dimensions d where each weak learner will correspond to a SDA. This randomization allows a single SDA weak learner to process pathologically correlated voxels across 3D local neighborhoods while still operating on the whole-brain image. Unlike the SVM objective which has a global optimum, SDAs can only converge to a local optimum via SG [1].We compensate for this by including a second level of randomization that samples sets of hyperparameters from a given hyper–parameter space. This basic structure drives the performance of our randomized denoising autoencoder (rDA).

Let 𝒱 = {1, ⋯, d} denote the indices of dimensions/voxels, and τ (υ), a distribution over υ ∈ 𝒱. In the simplest case, this can be a uniform distribution. We generate a bootstrap sample of B “blocks” where each block corresponds to input data along s_b dimensions/voxels (length of the block, fixed a priori). The mapping between voxels and blocks is given by τ (υ). Note that blocks may not be mutually exclusive (a voxel may belong to multiple blocks). Each block will be presented to T weak learners. Each of these weak learners correspond to a unique θ_t ∈ Θ for t = 1, …, T where Θ is the given hyper–parameter space. This means that each sample from the hyper-parameter space yields a weak learner. Our weak learner module is a L-layered stacked denoising autoencoder (SDA). The overall rDA architecture is an ensemble of B × T SDAs. Alg. 2.2 summarizes the block–wise training of rDA. Given training data as ${\{{\mathbf{x}}_i,{\mathbf{y}}_i\}}_1^n$, we first learn the transformations $({\mathbf{W}}_{b,t}^l,p_{b,t}^l,q_{b,t}^l)\forall b,t,l$. Denoting the L^th layer outputs by ${\mathbf{H}}_i={[[{\mathbf{z}}_{b,t}{\mathbf{h}}_{b,t}^L]]}_{1,1}^{B,T}$, the weighted regression pooling gives

where U are the regression coefficients and λ is the regularization constant. z is the known weight vector on B × T estimated L^th layer outputs. [[·]] denotes column–wise concatenation. The prediction for a new test input x is

The simplest choice for the block–wise sampler τ (υ) (i.e., at b = 1) assigns uniform probability over all dimensions/voxels. However, we can assign large weights on local neighborhoods which are more sensitive to dementia progression, if desired. Since d is large, we modify τ (υ) after each iteration (Reweigh step in Alg. 2.2) to prevent starvation of the previously unsampled dimensions. We can also setup τ (·) based on entropy or the result of a hypothesis test. Each weak learner output ${\mathbf{h}}_{b,t}^L$ is an estimate of y. Recall that SG learning is sensitive to the choice of hyper–parameters θ_t, particularly, the number of epochs and gradient learning rate influences the range and variance of these estimates [1,3]. Hence, randomization over θ_t mitigates this dependency by averaging over T such estimates for each block (i.e., set of dimensions/voxels). We can pool the block estimates via various means – average, using a ridge regression or other sophisticated schemes. But since SDAs are already capable of learning complex concepts [1], we use a simple linear combination with ℓ₂-loss providing minimum mean squared error. This addresses our goal of reducing the stochastic error of final predictions. Observe that rDA extends easily to multi–modal inputs by first constructing individual blocks for each modality and pooling across all modalities.

The sigmoid non–linearity ensures that rDA outputs are in ∈ [0, 1]. By labeling healthy controls as 1 and AD subjects as 0, we can then project the dementia scale to [0, 1]. The pooled outputs, referred to as rDA measure (rDAm), are then imaging–derived continuous predictors. We can then compute the sample sizes using rDAm as a customized outcome [11]. Denoting the mean change of rDA for placebo and treatment groups by μ_p and μ_t respectively, the number of subjects per arm is given by 2(Z_α/2+Z_1−β)²σ²/(μ_p−μ_t)² where 1−β is the desired power at significance level α. Using a conversion rate of ρ ∈ [0, 1] from MCI to AD, and inducing a drug effect of η (i.e., the treatment decreases the mean change by a fraction η), the sample size expression then simplifies to 2c²(Z_α/2+Z_1−β)²/(ηρ)² where c = σ/μ is the coefficient of variation. Since we only use one time–point data, the proportion ρ is set based on information from previously reported studies [13] Since rDA is an ensemble designed to reduce the prediction variance (and hence c), we hope to see much smaller sample sizes compared to others.

3.1 Data and Setup

We used Amyloid, FDG–PET and MRI data at baseline for 447 subjects (210 male, 237 female) from ADNI2 (Alzheimer’s Disease Neuroimaging Initiative). 131 were healthy (CN), 92 were demented (AD), 120 and 104 had early and late MCI (EMCI and LMCI) respectively. The labeling of EMCI and LMCI (done by ADNI) is based on the cognitive status of each subject. Of the 224 MCIs, 100 had maternal family history (FH) of AD, 52 had paternal and 23 had both. Pre-processing included extracting grey matter in normalized space, and correcting PET for average intensities in ponsvermis (FDG) and cerebellum (Amyloid). We train rDA on ADs (labeled 0) and CNs (labeled 1) alone, and test on MCIs. We use a multi-modal (MKL) ε–support vector regression (εMKm) as the baseline learning model [7]. Firstly, we evaluate if rDAm differentiates EMCI from LMCI. Additionally, we evaluated parental family history as a contributing risk factor. Since rDAm is a continuous marker, its correlations with CSF levels –τ, pτ, Aβ, τ/Aβ and pτ/Aβ (τ: τ-protein, pτ: phospho τ-protein, Aβ: Amyloid-β, are the cerebrospinal fluid protein levels, and very sensitive biomarkers of AD) are also assessed to verify if it is meaningful. We then estimate sample sizes using MCI to AD conversion rate of 37.7% [13]. rDA hyper-parameters in our experiments are L = 2, B = 1000 and T = 9, with uniform weighting (i.e., $z_{b,t}=\frac{1}{BT}\forall b,t$).

3.2 Results and Discussion

Table 1(a) and (b) show that rDAm is highly sensitive to EMCI vs. LMCI and the influence of FH. Although the baseline εMKm picks up these group differences, rDAm has much higher delineation power. In particular, the p–values for rDAm for FH positive vs. negative case are an order of magnitude smaller than that of εMKm. These show that rDAm is at least as good as a current state-of-the-art machine learning derived measures. Table 2 shows that rDAm has significant correlation (higher than εMKm in all but two cases) to CSF levels, which are proven biomarkers for MCI to AD progression [15]. Note that a negative correlation with say τ implies that rDAm decreases and the subject gets demented as τ increases. Specifically, higher correlations (with p ≤ 0.01) with pτ and pτ/Aβ suggest that rDAm is a useful continuous predictor. In most cases, these significance levels increase as more modalities are combined.

Table 3 shows that the coefficient of variation (CV) of rDAm for three different populations of interest – all MCIs, LMCIs and MCIs with positive FH. Observe that rDAm’s CV is smaller than that of εMKm for all three populations, and all possible combinations of modalities – making it a better candidate to be used as a prediction measure. Also, the CVs for MCIs with positive FH are smaller than that of late MCIs. This also suggests that a significant number of late MCIs currently have only a mild dementia in terms of rDAm. Fig. 1 shows the estimates on the same three populations of interest as above at 80% power (Refer to the supplement for more plots). Following Table 3, it should be straight forward to expect smaller sample estimates for rDAm compared to εMKm, which is exactly the case as shown in Fig. 1. Particularly, FDG and MRI gave smaller estimates than that of Amyloid following their smaller CV, reflecting that ~ 30% of healthy elderly have positive Amyloid scans. FH positive MCIs (last plot in Fig. 1) lead to much smaller sizes compared to using all MCIs and late MCIs. To get a sense of the improvement with respect to non–imaging based markers, we compared the best estimate (over all modalities) of εMKm and rDAm with that of MMSE and CSF levels in Table 4. At 80% power, the best estimates across CSF markers was 973 and 975 (for τ and τ/Aβ respectively) compared to that of 193 using rDAm – more than 5–fold decrease. It should be noted that all these estimates use only “single time–point” data combined with known conversion rates, in contrast to direct longitudinal measurement [8,4]. Hence, the sizes using MMSE and CSF are as high as 1500, indicating that estimates on the order of two hundred (that of rDAm) are highly significant. Overall, the results show that imaging–derived markers lead to much smaller trials than cognitive scores and/or CSF levels.