We synthesize multivariate correlation and regression methods to characterize unique relationships among compositional and physical properties of a planetary surface locally, regionally, and globally. Martian data including elemental mass fractions, areal fractions of mineral types, and thermal inertia constitute our case study. We incorporate techniques to address the effects of spatial autocorrelation and heteroscedasticity. We also utilize method and fit diagnostics. While the Mars Odyssey and Mars Global Surveyor missions provide the exploratory context in our discussion, our approach is applicable whenever the interrelationships of spatially binned data of continuous-valued planetary attributes are sought. For example, our regional-scale case study reinforces the strength of the spatial correlation among K, Th, and the dominant mineralogic type within northern low albedo regions (surface type 2) of Mars. Recent chemical and mineralogic data from the MESSENGER mission at Mercury and Dawn at Vesta may be analyzed effectively with these hierarchical regression methods to constrain geochemical processes. Likewise, our algorithm could be applied locally with the wide variety of compositional data expected from the MSL mission at Gale Crater in general, and the ChemCam sampling grids in particular.