Regression analysis for areal data is common in numerous fields, including public health, ecology, and econometrics. Often, the goal of such an analysis is to quantify relationships between an outcome or outcomes of interest and covariates. In our present work, we propose several approaches to modeling areal data in areas including neuroimaging, cancer epidemiology, and demography. Much of our work is driven by the need to efficiently model large datasets with spatial dependencies. For instance, we model functional magnetic resonance imaging (fMRI) data using spatial Bayesian variable selection for detecting task-based brain activity. Fitting a full statistical model to large fMRI datasets can be computationally burdensome, so we efficiently implement a full statistical model via a brain parcellation that allows for parallel computations. Brain activity within each parcel is modeled as regressions located on a lattice with latent indicators for regressors, representing zero or non-zero activity levels. A sparse areal mixed model (SAMM) captures spatial dependence among indicator variables for a given stimulus. The SAMM permits more efficient computation than traditional spatial models typically employed. Through simulation we show that our parcellation scheme performs well in various scenarios. We apply our method to data from a task-based fMRI experiment. Large multivariate and zero-inflated spatial data motivates another area where we develop a computationally efficient model. Our proposed multivariate model employs the SAMM to effectively reduce the dimensionality of the spatial effects and alleviate bias and variance inflation that affects estimation involving fixed effects covariates. We then adapt the multivariate model to handle zero-inflated count data via the truncated Poisson hurdle model. We carry out several simulation studies to evaluate the importance of modeling the correlation between the spatial effects of the different outcomes. We apply our new method on a high-dimensional zero-inflated dataset collected as part of the 2010 US Census. Another motivation for our work is the conditional autoregressive (CAR) model and its numerous drawbacks. We develop a covariance selection model that has a marginal interpretation that eases the restriction of interpreting covariate effects conditionally. Additionally, though the proper CAR model's dependence parameter has an intuitive conditional interpretation, the marginal interpretation is complicated and counterintuitive. To overcome these drawbacks, we introduce a copula-based covariance selection model. We achieve unbiased estimation of marginal parameters with an intuitive marginal interpretation. The covariance selection copula's single dependence parameter is the first-order correlation which provides a dependence structure with intuitive marginal interpretations. We develop a computational framework that permits efficient frequentist inference. We evaluate small and large sample performance of our method under simulated conditions and apply our procedure to the widely studied Slovenia stomach cancer dataset.