A Bayesian Convolutional Neural Network-based Generalized Linear Model

A Bayesian perspective can quantify uncertainties through posterior densities of $\bm{\theta}$. Specifically, Gal and Ghahramani (2016a) proposed an efficient Bayesian alternative by regarding weight parameters in CNN’s kernels as random variables. By extending a result in Gal and Ghahramani (2016b), Gal and Ghahramani (2016a) showed that applying dropout after every convolution layer can approximate the intractable posterior distribution of deep GP.

Let $\mathbf{X}_{n}\in{\mathbb{R}}^{M_{0}\times R_{0}\times C_{0}}$ be an input with $M_{0}$ height, $R_{0}$ width, and $C_{0}$ channels. For example, $M_{0}$ and $R_{0}$ are determined based on the resolution of images, and $C_{0}=3$ for color images ($C_{0}=1$ for black and white images). There are three types of layers in CNN: the convolution, pooling, and fully connected layers. Convolution layers move kernels (filters) on the input image and extract information, resulting in a “feature map”. Each cell in the feature map is obtained from element-wise multiplication between the image and the filter matrix. The weight values in filter matrices are parameters to be estimated. Higher output values indicate stronger signals on the image, which can account for spatial structures. We can extract the important neighborhood features from the input by shifting the kernels over all pixel locations with a certain step size (stride). Then, the pooling layers are applied to downsample the feature maps. The convolution and pooling layers are applied repeatedly to reduce the dimension of an input image. Then, the extracted features from the convolution/pooling layers are flattened and connected to the output layer through a fully connected layer, which can learn non-linear relationships between the features and the output.

Consider a CNN with $L$ number of layers, where the last two are fully connected and output layers. For the $l$th convolution layer, there are $C_{l}$ number of kernels $\mathbf{K}_{l,c}\in{\mathbb{R}}^{H_{l}\times D_{l}\times C_{l-1}}$ with $H_{l}$ height, $D_{l}$ width, and $C_{l-1}$ channels for $c=1,\cdots,C_{l}$. Note that the dimension of the channel is determined by the number of kernels in the previous convolution layer. After applying the $l$th convolution layer, the $(m,r)$th element of the feature matrix $\bm{\eta}_{n,(l,c)}\in{\mathbb{R}}^{M_{l}\times R_{l}}$ becomes

where $b_{l,c}$ is the bias parameter. For the first convolution layer, $\bm{\eta}_{n,(l-1,c)}$ is replaced with an input image $\mathbf{X}_{n}$. (i.e., $\bm{\eta}_{n,(0,c)}=\mathbf{X}_{n}$). For computational efficiency and stability, ReLU functions are widely used as activation functions $\sigma(\cdot)$. By shifting kernels, the convolution layer can capture neighborhood information, and the neighborhood structure is determined by the size of the kernels. Gal and Ghahramani (2016a) pointed out that extracting the features from the convolution operation can be reformulated as an affine transformation in standard DNN (details are in Web Appendix A.2). Therefore, CNNs can also be represented as a deep GP with dropout approximation; applying dropout after every convolution layer before pooling can approximate a deep GP.

Note that optimal hyperparameter tuning in NNs requires a number of empirical fittings by varying hyperparameters. Wang et al. (2019) pointed out that a higher dropout rate can lead to a slower convergence rate, while a lower rate can deteriorate generalization performance. Although Gal and Ghahramani (2016b, a) used a dropout rate of around 0.1-0.5, it does not always guarantee optimal performance. Therefore, we recommend conducting empirical studies by varying dropout rates within such a range for validation. In Web Appendix B, we observe that our method is robust across different dropout rates.

After the $L-2$th convolution layers, the extracted feature tensor $\bm{\eta}_{n,L-2}\in{\mathbb{R}}^{M_{L-2}\times R_{L-2}\times C_{L-2}}$ is flattened in the fully connected layer. By vectorizing $\bm{\eta}_{n,L-2}$, we can represent the feature tensor as a feature vector $\bm{\phi}_{n,L-2}=\operatorname{vec}(\bm{\eta}_{n,L-2})\in{\mathbb{R}}^{k_{L-2}}$, where $k_{L-2}=M_{L-2}R_{L-2}C_{L-2}$. Finally, we can extract a feature vector $\bm{\phi}_{n,L-1}\in{\mathbb{R}}^{k_{L-1}}$ from the fully connected layer. From here, we define the collection of the extracted feature vector as $\bm{\Phi}=\{\bm{\phi}_{n,L-1}\}_{n=1}^{N}$. We suppress the index $L-1$ in $\bm{\Phi}_{L-1}$ for notational convenience.

The training step of BayesCNN is identical to that of standard CNN. In the prediction step, the algorithm performs forward propagation for the given estimated NN parameters. For each iteration, $\{\mathbf{d}_{l}\}_{l=2}^{L}$ are randomly sampled from the Bernoulli distribution with rate $\psi_{l}$ and are applied to the features. The algorithm repeats this $M$ times to obtain an MC sample of $M$ predicted outputs. We provide implementation details for BayesCNN in Web Appendix C.1.

We propose BayesCGLM by regressing $\mathbf{Y}$ on $\bm{\Phi}$ and additional covariates $\mathbf{Z}$. Here, we use $\bm{\Phi}$ as a basis matrix that encapsulates information of high-dimensional and correlated input variables $\mathbf{X}$. Let $\mathbf{A}=[\mathbf{Z},\bm{\Phi}]\in{\mathbb{R}}^{N\times(p+k_{L-1})}$ be the design matrix. Then, our model is

where $\bm{\beta}=(\bm{\gamma}^{\top},\bm{\delta}^{\top})^{\top}\in{\mathbb{R}}^{p+k_{L-1}}$ is the corresponding regression coefficients and $g(\cdot)$ is a one-to-one continuously differential link function. The proposed method can provide an interpretation of the regression coefficients and quantify uncertainties in prediction and estimation. Furthermore, the training step does not require additional computational costs compared to the traditional CNN. We train BayesCNN with the negative log-likelihood loss function corresponding to the type of output variable. To be more specific, we used the mean squared loss, binary cross-entropy loss, and Poisson loss for Gaussian, binary, and count outputs, respectively. Such choices are natural in that the loss functions are directly defined by the distribution of the output variable.

Since we place a posterior distribution over the kernels of the fitted BayesCNN (Section 3.1), we can generate MC samples of the features $\bm{\Phi}^{(m)}$, $m=1,\cdots,M$ from the variational distribution. To account for the uncertainties in estimating the features, we use all the MC samples of the features rather than their point estimates. For each MC sample of $\mathbf{A}^{(m)}=[\mathbf{Z},\bm{\Phi}^{(m)}]$, we can fit individual GLM in parallel; we have $M$ number of posterior distributions of $\bm{\beta}_{m}=(\bm{\gamma}_{m}^{\top},\bm{\delta}_{m}^{\top})^{\top}\in{\mathbb{R}}^{p+k_{L-1}}$. This step is parallelizable, so computational wall times decrease with more cores. Then, we construct an aggregated distribution (the “ensemble-posterior”) of $\bm{\beta}$. In the following section, we provide details about constructing the ensemble-posterior distribution. The fitting procedure of BayesCGLM is summarized in Web Appendix C.2.

In general, inference for the regression coefficients $\bm{\gamma}$ in (8) is challenging for DNN. Although one can regard the weight parameters at the last hidden layer as $\bm{\gamma}$, we cannot obtain uncertainties in estimates. Even BayesDNN (Gal and Ghahramani, 2016b) only provided point estimates of weight parameters, while they can quantify uncertainties in response prediction. On the other hand, our method directly provides the posterior distribution of $\bm{\gamma}$, which is useful for standard inference. In addition, the complex structure of DNN often leads to nonconvex optimization (Goodfellow et al., 2014; Dauphin et al., 2014), resulting in an inaccurate estimate of $\bm{\gamma}$. Since the stopping rule of the SGD algorithm is based on prediction accuracy (not on convergence in parameter estimation), the convergence of individual parameters, including $\bm{\gamma}$ cannot be guaranteed. The problem gets further complicated by a large number of parameters that need to be simultaneously updated in the SGD algorithm. Note that BayesCGLM also uses the SGD algorithm to extract $\bm{\Phi}$. From this, we can obtain one of the optimal combinations of weight parameters and the corresponding $\bm{\Phi}$. Although such $\bm{\Phi}$ may not be the global optimizer, it is still informative to predict responses because it is obtained by minimizing the loss function. However, this does not guarantee the convergence of individual weight parameters in nonconvex optimization. Instead, our two-stage approach can alleviate this issue by using extracted $\bm{\Phi}$ as a fixed basis matrix in GLM. In Web Appendix D, we conduct the simulation studies under different model configurations and show that BayesCGLM can recover the true $\bm{\gamma}$ values well, while estimates from BayesCNN can be biased.

Here, we describe a Bayesian inference for the proposed method. We first define the posterior distribution of the regression coefficient $\pi(\bm{\beta}|\mathbf{D})$. Let $\{\mathbf{W}_{l},\mathbf{b}_{l}\}_{l=1}^{L-1}$ be the weights and biases corresponding to all hidden layers in CNN. We can vectorize kernels in CNN and represent them as weights as in the standard DNN (see Web Appendix A.2 for details). Note that $\bm{\Phi}$ can be deterministically obtained from forward propagation for the given $\mathbf{X}$ and $\{\mathbf{W}_{l},\mathbf{b}_{l}\}_{l=1}^{L-1}$. Therefore, the posterior distribution can be represented as

where, $\pi(\bm{\beta}|\mathbf{D},\{\mathbf{W}_{l},\mathbf{b}_{l}\}_{l=1}^{L-1})$ is the conditional posterior, and $\pi(\mathbf{W}_{l}|\mathbf{D})$, $\pi(\mathbf{b}_{l}|\mathbf{D})$ are marginal posteriors for weight and bias, respectively. Since it is challenging to compute (9) directly, we approximate it through MC dropout as

where $q(\mathbf{W}_{l})$ and $q(\mathbf{b}_{l})$ are variational distributions defined in (4). Then the MC approximation to (10) is

Here $\{\{\mathbf{W}^{(m)}_{l},\mathbf{b}^{(m)}_{l}\}_{l=1}^{L-1}\}_{m=1}^{M}$ are sampled from (4). One may choose the standard MCMC algorithm to estimate each $\pi(\bm{\beta}_{m}|\mathbf{D},\{\mathbf{W}^{(m)}_{l},\mathbf{b}^{(m)}_{l}\}_{l=1}^{L-1})$ in (11). However, this may require a long run of the chain to guarantee good mixing, resulting in excessive computational cost. Instead, we use the Laplace method to approximate each $\pi(\bm{\beta}_{m}|\mathbf{D},\{\mathbf{W}^{(m)}_{l},\mathbf{b}^{(m)}_{l}\}_{l=1}^{L-1})$. In our preliminary studies, we observed that the Laplace approximation leads to a similar inference result as the MCMC method within a much shorter computing time when the sample size is large enough.

For the Laplace approximation, we first compute $\bm{\Phi}^{(m)}$ for the given $\{\mathbf{W}^{(m)}_{l},\mathbf{b}^{(m)}_{l}\}_{l=1}^{L-1}$ through forward propagation. Then, we obtain the maximum likelihood estimate (MLE) $\widehat{\bm{\beta}}_{m}$ using GLM by regressing $\mathbf{Y}$ on $\bm{\Phi}^{(m)}$ and $\mathbf{Z}$. We approximate the posterior of $\bm{\beta}_{m}$ as $\mathcal{N}(\widehat{\bm{\beta}}_{m},\widehat{\mathbf{B}}^{-1}_{m})$, where $\widehat{\mathbf{B}}_{m}\in{\mathbb{R}}^{(p+k_{L-1})\times(p+k_{L-1})}$ is the observed Fisher information matrix from the $m$th MC samples. From this procedure, (11) can be approximated as

where $\varphi(\bm{x};\bm{\mu},\bm{\Sigma})$ is a multivariate normal density with mean $\bm{\mu}$ and covariance $\bm{\Sigma}$. In summary, (12) is the approximation of the posterior distribution $\pi(\bm{\beta}|\mathbf{D})$ through MC dropout and the Laplace method. Note that BayesCGLMs approximate $\pi(\bm{\beta}|\mathbf{D})$ through the MC average from the same probability model (i.e., data generating mechanism). Therefore, the interpretation of $\bm{\beta}$ does not change from ensemble to ensemble. However, the interpretation of $\bm{\delta}$ can be more difficult because each column of $\bm{\Phi}$ comes from complex (and repeated) convolution operations. It is challenging to explain the impact of changes in features based on the estimated $\bm{\delta}$.

Similar to BayesCNN, our method can quantify uncertainties in predictions, which is an advantage over deterministic NNs. Note that BayesCGLM can also provide interpretations of model parameters with uncertainty quantification, while BayesCNN cannot. For unobserved response $\mathbf{Y}^{\ast}\in{\mathbb{R}}^{N_{\text{test}}}$, we have $\mathbf{A}^{\ast}_{m}=[\mathbf{Z}^{\ast},\bm{\Phi}^{\ast(m)}]\in{\mathbb{R}}^{N_{\text{test}}\times(p+k_{L-1})}$ for $m=1,\cdots,M$. Here, $\bm{\Phi}^{\ast(m)}$ is $m$th features extracted from the last hidden layer of the trained BayesCNN for given $\mathbf{X}^{\ast}$ and $\mathbf{Z}^{\ast}$. Then, from (12), the predictive distribution of the linear predictor is

For given (13), we can generate $\mathbf{Y}^{\ast}$ from the underlying distribution. For example, when we have a Gaussian response, we have $\mathbf{Y}^{\ast}\sim\mathcal{N}(\frac{1}{M}\sum_{m=1}^{M}\mathbf{A}^{\ast}\bm{\widehat{\beta}}_{m},\widehat{\sigma}^{2})$ where $\widehat{\sigma}^{2}=\sum_{m=1}^{M}(\mathbf{A}^{(m)}\bm{\widehat{\beta}}_{m}-\mathbf{Y})^{\top}(\mathbf{A}^{(m)}\bm{\widehat{\beta}}_{m}-\mathbf{Y})/NM$. For the count response, we can simulate $\mathbf{Y}^{\ast}$ from the Poisson distribution with an intensity of $\exp(\frac{1}{M}\sum_{m=1}^{M}\mathbf{A}^{\ast}\bm{\widehat{\beta}}_{m})$.

As we described, BayesCGLM involves several approximation steps: (1) MC dropout and (2) the Laplace method; therefore, investigating approximation quality is crucial. In Web Appendix D.5, we conduct a numerical study under the simple NN setting. Here, we use a standard NN to extract features; therefore, we refer to our method as BayesDGLM. We compare the performance of BayesDGLM, BayesDNN, and a fully Bayesian method (MCMC). We observe that the BayesDGLM provides comparable inference results to the fully Bayesian method with much smaller computational costs. Furthermore, BayesDNN cannot quantify uncertainties in parameter estimates, while BayesDGLM can. Considering that it is infeasible to apply the fully Bayesian method for complex DNNs, BayesDGLM would be a practical option in many applications. We provide details for the experiment in Web Appendix D.5.

According to Noone et al. (2020), brain and nervous system tumors have low 5-year relative survival rates, which have gradually increased from 1975 to 2016, reaching 33%. Furthermore, some brain tumors can also be cancerous. Considering that the 5-year survival rate for overall cancers is 67%, brain tumors can be life-threatening. MRI is the standard diagnostic tool for brain tumors, and texture-based first- and second-order statistics of MRI images are useful for classification (Aggarwal and Agrawal, 2012), and these statistics are widely used in the classification of biomedical images (Ismael and Abdel-Qader, 2018). However, these statistics do not capture the spatial information of correlated structures in MRI images. Automated methods for brain tumor classification are becoming increasingly important, as manual analysis of MRI data by medical experts is challenging. It is also crucial to convey the uncertainty associated with these automated methods to ensure statistically sound decision-making. Compared to previous studies (Aggarwal and Agrawal, 2012; Ismael and Abdel-Qader, 2018), we use both MRI images and texture-based statistics to predict brain tumors while quantifying uncertainties. The dataset is collected from the Brain Tumor Image Segmentation Challenge (Menze et al., 2014). The binary response indicates whether a patient has a fatal brain tumor. For each patient, we have an MRI image with its two summary variables (first and second-order features of MRI images). Here, we use $N=2,508$ observations to fit our model and reserve $N_{\text{test}}=1,254$ observations for validation. We fit BayesCNN using $240\times 240$ pixel gray images $\mathbf{X}$, scalar covariates $\mathbf{Z}$ as input, and the binary response $\mathbf{Y}$ as output. From the fitted BayesCNN, we can extract $\bm{\Phi}\in{\mathbb{R}}^{2,508\times 16}$ from the last hidden layer. In this study, we use two summary variables (first and second-order features of MRI images) as covariates. With the extracted features from the BayesCNN, we fit logistic GLMs by regressing $\mathbf{Y}$ on $[\mathbf{Z},\bm{\Phi}^{(m)}]\in{\mathbb{R}}^{2,508\times 18}$ for $m=1,\cdots,500$.

Table 1 shows inference results from different methods. Our method shows that the first-order feature has a negative relationship, while the second-order feature has a positive relationship with the brain tumor risk. Both coefficients are statistically significant based on the highest posterior density (HPD) intervals. Note that previous work (Aggarwal and Agrawal, 2012) only predicted tumor status but did not provide such an interpretation from the first and second-order statistics. Therefore, we can conclude that darker MRI images (larger first-order statistics) with lower variation (smaller second-order statistics) indicate a lower possibility of tumor presence. We observe that the sign of $\bm{\gamma}$ estimates from BayesCGLM is aligned with those from GLM.

Figure 1 is a score plot with the first and second principal components of $\bm{\Phi}$. Figure 1 (a) shows a clear separation of tumor and non-tumor status based on two principal components. In particular, we observe that patients without tumors mostly have positive values of the first principal component (62.9% explained variability). This indicates that $\bm{\Phi}$ contains useful information from images for predicting brain tumor status. By incorporating $\bm{\Phi}$ information, BayesCGLM shows the most accurate prediction performance while providing credible intervals for the estimates.

Figure 2 (a) shows that the area under the curve (AUC) for the receiver operating characteristic curve (ROC) plot is about 0.96, which shows outstanding prediction performance. The predicted probability surface becomes larger for the binary responses with a value of 1 (Figure 2 (b)). Furthermore, the true and predicted binary responses are well aligned. We also investigate how prediction uncertainties are related to actual images. The top panel in Figure 2 (c) shows four examples of correctly classified cases with low prediction uncertainties. The results show that when the image is overall dark, and the tumor area clearly stands out from the background, we have low uncertainties. On the other hand, the bottom panel in Figure 2 (c) shows four examples of misclassified cases with high prediction uncertainties. In these cases, the background image is too bright, so the tumor area is hardly distinguishable from the background. Although the point predictions misclassify these cases, higher uncertainty informs us that further investigation is needed. Therefore, better medical decisions can be made when uncertainty measures are considered.

Anxiety disorders are becoming a prevalent mental health issue worldwide, with studies showing a more than 25% increase in cases in 2020 (Organization et al., 2022). There have been several works to study the relationship between anxiety scores and psychological assessments collected from the survey data (Kaplan et al., 2015). However, no previous studies considered the fMRI data, which contains important information about the functional connectivity of the brain. Here, we utilize fMRI information through the extracted feature to construct an anxiety score model. From this, we can incorporate the correlation structure of the brain when we study the relationships between psychological assessments and anxiety levels.

The resting-state fMRI data is collected by the enhanced Nathan Kline Institute-Rockland Sample (NKI-RS) to create a large-scale community sample of participants across the lifespan. For 156 patients, resting-state fMRI signals were scanned for 394 seconds. Then we transform fMRI images into 278 functional brain Regions-of-Interest (ROIs) following the parcellation scheme of Shen et al. (2013).

A functional network was computed using the pairwise Pearson correlation coefficient, constructing $278\times 278$ square symmetric matrix, which is called the ROI-signal matrix. The correlation matrix of brain functions has been widely used to study the functional connectivity of the brain (Sporns, 2002). Therefore, we also use the correlation matrix as input $\mathbf{X}$. Participants in the study completed a survey that included various physiological and psychological assessments. In this study, we use $N=156$ observations to fit our model and reserve $N_{\text{test}}=52$ observations for validation. We use conscientiousness, extraversion, agreeableness, and neuroticism as $\mathbf{Z}$ and the averaged anxiety score as $\mathbf{Y}$. We train BayesCNN with $\mathbf{X}$, $\mathbf{Z}$ as input, and the continuous response $\mathbf{Y}$ as output. From the trained BayesCNN, we extract each $m$th feature matrix $\bm{\Phi}^{(m)}\in{\mathbb{R}}^{104\times 8}$ from the last hidden layer. Then we regress $\mathbf{Y}$ on $[\mathbf{Z},\bm{\Phi}^{(m)}]$ for $m=1,\cdots,500$.

Table 2 indicates that BayesCGLM has comparable prediction performance compared to BayesCNN and can quantify uncertainties of regression coefficients. While the sign of the estimated $\bm{\gamma}$ are identical across different methods, there are differences in the statistical significance, particularly in agreeableness and conscientiousness variables. From BayesCGLM, we can infer that the person who has emotional instability and vulnerability to unpleasant emotions (neuroticism) and who is less likely to rebel against others (agreeableness) has a high anxiety and depression score, which are aligned with the previous study in Kaplan et al. (2015). On the other hand, individuals who are more sociable (extroverted) and driven to achieve their goals (conscientiousness) tend to exhibit lower anxiety and depression scores, as also supported by the previous studies (Naragon-Gainey et al., 2014; Fan et al., 2020). Note that these previous studies neither studied the psychological assessments simultaneously nor fMRI information.

GLM, which does not utilize fMRI information, draws different conclusions in a significance test. Table 2 shows that agreeableness and conscientiousness variables are not statistically significant in GLM, which contradicts the findings in the previous studies (Kaplan et al., 2015; Fan et al., 2020). Figure 1 is a score plot with the first and second principal components of $\bm{\Phi}$. Figure 1 (b) indicates a negative correlation between anxiety scores and the first principal component (62.6% explained variability). Therefore, $\bm{\Phi}$ encapsulates information about brain functions that are highly associated with anxiety levels. By incorporating such important information into the model, BayesCGLM can accurately infer the relationship between psychological assessments and show improved prediction performance.

We observe that all the methods provide similar coefficient estimates. We can infer that the person who has emotional instability and vulnerability to unpleasant emotions (neuroticism) and who is less likely to rebel against others (agreeableness) has a high anxiety and depression score, which gives the comparable result from Kaplan et al. (2015). The person who is more inclined to seek sociability and is eager to attain goals (extraversion), on the other hand, has lower anxiety and depression scores. Note that the previous work (Kaplan et al., 2015) only conducted exploratory data analysis based on correlation coefficients without quantifying the uncertainty. On the other hand, our model can estimate the coefficients by incorporating correlated fMRI information.

Figure 3 shows that there is some agreement between the true and predicted responses. We observe that HPD prediction intervals include the true response well (92.3% prediction coverage). Due to the small sample size, prediction intervals are wider than previous examples, which is natural. In Figure 3, triangles represent individuals with the three largest HPD intervals. We observe that the anxiety score is underestimated for these people, but the high uncertainty informs us that we may need to pay additional attention to these individuals with severe anxiety issues; they might have been overlooked without uncertainty quantification.