Medical imaging technologies have been generating extremely complex data sets. This dissertation makes further contributions to the development of statistical tools motivated by modern biomedical challenges. Specifically we develop methods to characterize varying associations between ultra-high dimensional imaging data and low-dimensional clinical outcomes. The first part of this dissertation is motivated by the major limitations faced by traditional voxel-wise models, where voxels are commonly treated as independent units, and the assumption of Gaussian distribution of the neuroimaging measurements is usually flawed. We develop a class of hierarchical spatial transformation models to model the spatially varying associations between imaging measurements in a three-dimensional (3D) volume (or 2D surface) and a set of covariates. The proposed approach include a spatially varying Box-Cox transformation model and a Gaussian Markov random field model. The second part is motivated by the challenges faced by ultra-high dimensional datasets. In particular, we introduce a method to predict clinical outcomes from ultra-high dimensional covariates. The proposed models reduce dimensionality to a manageable level and further apply dimension reduction techniques, e.g. principal components analysis and tensor decompositions to extract and select low-dimensional important features.