We develop a novel statistical model, called multiscale adaptive regression model (MARM), for spatial and adaptive analysis of neuroimaging data. The primary motivation and application of the proposed methodology is statistical analysis of imaging data on the two-dimensional (2D) surface or in the 3D volume for various neuroimaging studies. The existing voxel-wise approach has several major limitations for the analyses of imaging data, underscoring the great need for methodological development. The voxel-wise approach essentially treats all voxels as independent units, whereas neuroimaging data are spatially correlated in nature and spatially contiguous regions of activation with rather sharp edges are usually expected. The initial smoothing step before the voxel-wise approach often blurs the image data near the edges of activated regions and thus it can dramatically increase the numbers of false positives and false negatives. The MARM, which is developed for addressing these limitations, has three key features in the analysis of imaging data: being spatial, being hierarchical, and being adaptive. The MARM builds a small sphere at each location (called voxel) and use these consecutively connected spheres across all voxels to capture spatial dependence among imaging observations. Then, the MARM builds hierarchically nested spheres by increasing the radius of a spherical neighborhood around each voxel and combine all the data in a given radius of each voxel with appropriate weights to adaptively calculate parameter estimates and test statistics. Theoretically, we first establish that the MARM outperforms classical voxel-wise approach. Simulation studies are used to demonstrate the methodology and examine the finite sample performance of the MARM. We apply our methods to the detection of spatial patterns of brain atrophy in a neuroimaging study of Alzheimer's disease. Our simulation studies with known ground truth confirm that the MARM significantly outperforms the voxel-wise methods.