To study the evolution of several quantitative traits, the classical phylogenetic comparative framework consists of a multivariate random process running along the branches of a phylogenetic tree. The Ornstein-Uhlenbeck (OU) process is sometimes preferred to the simple Brownian motion (BM) as it models stabilizing selection toward an optimum. The optimum for each trait is likely to be changing over the long periods of time spanned by large modern phylogenies. Our goal is to automatically detect the position of these shifts on a phylogenetic tree, while accounting for correlations between traits, which might exist because of structural or evolutionary constraints. We show that, in the presence of shifts, phylogenetic Principal Component Analysis fails to decorrelate traits efficiently, so that any method aiming at finding shifts needs to deal with correlation simultaneously. We introduce here a simplification of the full multivariate OU model, named scalar OU, which allows for noncausal correlations and is still computationally tractable. We extend the equivalence between the OU and a BM on a rescaled tree to our multivariate framework. We describe an Expectation-Maximization (EM) algorithm that allows for a maximum likelihood estimation of the shift positions, associated with a new model selection criterion, accounting for the identifiability issues for the shift localization on the tree. The method, freely available as an R-package (PhylogeneticEM) is fast, and can deal with missing values. We demonstrate its efficiency and accuracy compared to another state-of-the-art method ($\ell$1ou) on a wide range of simulated scenarios and use this new framework to reanalyze recently gathered data sets on New World Monkeys and Anolis lizards.