# A Novel Bayesian Linear Regression Model for the Analysis of Neuroimaging Data
## 3 Methods
### A Dual Bayesian Linear Regression Model with Feature Selection
#### 3.1.2 Generative Model  
As the graphical model of Figure 2 shows, the generative model of DBL-FS considers that each datum, \(\textbf{x}_{n,i}\), is combined with a weight vector **w** plus some Gaussian noise to generate the output variable:  
\[y_{n}=\textbf{x}_{n,i}\textbf{w}+\eta, \tag{1}\]  
where \(\eta\) is a Gaussian noise with zero mean and precision \(\tau\). In turn, the noise precision is modeled with a gamma distribution with parameters \(a_{0}^{\tau}\), \(b_{0}^{\tau}\):  
\[\tau\sim\Gamma(a_{0}^{\tau},b_{0}^{\tau}) \tag{2}\]  
In addition, DBL-FS considers that the weight associated to the \(d\)-th input feature follows a normal distribution:  
\[w_{d}\sim\mathcal{N}\Big{(}0,\;\alpha_{d}^{-1}\Big{)}\quad d=1,\ldots,D \tag{3}\]  
where its precision, \(\alpha_{d}\), is modeled with a gamma distribution as:  
\[\alpha_{d}\sim\Gamma(a_{0}^{\alpha},b_{0}^{\alpha})\quad d=1,\ldots,D \tag{4}\]  
This ARD prior over \(w_{d}\) allows us to obtain the relevance over the elements of **w**, and therefore, DBL-FS is capable of automatically setting to zero the features that are irrelevant for the problem.  
As the model will have to work with MRI data, composed by few samples (less than 100) and tens or hundreds of thousands of voxels, it is clear that working in the primalspace is not the most efficient way to proceed. So, we propose to reformulate the model making use of the Representer Theorem [33] (RT). That is, as the RT states that the primal weights of any regression model resulting from minimizing an empirical error (risk) can be expressed as a linear combination of the input data and its equivalent dual coefficients, we can express \(\mathbf{w}\) as:  
\[\mathbf{w}=\mathbf{X}^{T}\mathbf{a} \tag{5}\]  
where \(\mathbf{a}\) is a vector of length \(N\) containing the dual variables. As we will see later (see Equation (19)), the lower bound that maximizes our variational inference is equivalent to minimizing an empirical cost. This way, the model can be formulated to work in the dual space as:  
\[y_{n}=\mathbf{k}_{n;\mathbf{a}}+\eta, \tag{6}\]  
where \(\mathbf{k}_{n;\mathbf{a}}\) denotes to the \(n\)-th row of the linear kernel matrix \(\mathbf{K}=\mathbf{X}\mathbf{X}^{T}\). This way, with the dual formulation, the target variables \(y_{n}\), for \(n=1,\ldots,N\), are modeled as:  
\[y_{n}\sim\mathcal{N}\Big{(}\mathbf{k}_{n;\mathbf{a}},\tau^{-1}\Big{)}\quad n= 1,\ldots,N. \tag{7}\]  
With this new formulation, the model will be able to work in the space of \(\mathbf{a}\), where only \(N\) parameters have to be inferred. Thus, model complexity and overfitting risks are drastically reduced, as long as we are able to maintain the feature relevance determination over \(\mathbf{w}\), providing the model with feature selection capabilities.  
Finally, it is important to note that the model formulation does not need to specifically include the distribution of \(\mathbf{a}\) since the relation between \(\mathbf{w}\) and \(\mathbf{a}\) is deterministic, and therefore, the statistical characterization of \(\mathbf{w}\) is also characterizing \(\mathbf{a}\).