  Task functional magnetic resonance imaging (fMRI) is a type of neuroimaging
data used to identify areas of the brain that activate during specific tasks or
stimuli. These data are conventionally modeled using a massive univariate
approach across all data locations, which ignores spatial dependence at the
cost of model power. We previously developed and validated a spatial Bayesian
model leveraging dependencies along the cortical surface of the brain in order
to improve accuracy and power. This model utilizes stochastic partial
differential equation spatial priors with sparse precision matrices to allow
for appropriate modeling of spatially-dependent activations seen in the
neuroimaging literature, resulting in substantial increases in model power. Our
original implementation relies on the computational efficiencies of the
integrated nested Laplace approximation (INLA) to overcome the computational
challenges of analyzing high-dimensional fMRI data while avoiding issues
associated with variational Bayes implementations. However, this requires
significant memory resources, extra software, and software licenses to run. In
this article, we develop an exact Bayesian analysis method for the general
linear model, employing an efficient expectation-maximization algorithm to find
maximum a posteriori estimates of task-based regressors on cortical surface
fMRI data. Through an extensive simulation study of cortical surface-based fMRI
data, we compare our proposed method to the existing INLA implementation, as
well as a conventional massive univariate approach employing ad-hoc spatial
smoothing. We also apply the method to task fMRI data from the Human Connectome
Project and show that our proposed implementation produces similar results to
the validated INLA implementation. Both the INLA and EM-based implementations
are available through our open-source BayesfMRI R package.
