In two articles, we present ‘coregionalization analysis with a drift’ (CRAD), a method to assess the multi-scale variability of and relationships between ecological variables from a multivariate spatial data set. In phase I of CRAD (the first article), a deterministic drift component representing the large-scale pattern and a random component modeled as a second-order stationary process are estimated for each variable separately. In phase II (this article), a linear model of coregionalization (LMC) is fitted by estimated generalized least squares to the direct and cross experimental variograms of residuals (i.e., after the removal of estimated drifts). Structural correlations and coefficients of determination at smaller scales are then computed from the estimated coregionalization matrices, while the estimated drifts are used to calculate pseudo coefficients at large scale. The performance of five procedures in estimating correlations and coefficients of determination was compared using a Monte Carlo study. In four CRAD procedures, drift estimation was based on local polynomials of order 0, 1, 2 (L0, L1, L2) or a global polynomial with forward selection of the basis functions; the fifth procedure was coregionalization analysis (CRA), in which large-scale patterns were modeled as a supplemental component in the LMC. In bivariate and multivariate analyses, the uncertainty in the estimation of correlations and coefficients of determination could be related to the interference between spatial components within a bounded sampling domain. In the bivariate case, most procedures provided acceptable estimates of correlations. In regionalized redundancy analysis, uncertainty was highest for CRA, while L1 provided the best results overall. In a forest ecology example, the identification of scale-specific correlations between plant species diversity and soil and topographical variables illustrated the potential of CRAD to provide unique insight into the functioning of complex ecosystems.