In many applications, interest focuses on assessing the relationship between a predictor and a multivariate outcome variable, and there may be prior knowledge about the shape of the regression curves. For example, regression functions that relate dose of a possible risk factor to different adverse outcomes can often be assumed to be nondecreasing. In such cases, interest focuses on (1) assessing evidence of an overall adverse effect, (2) determining which outcomes are most affected, and (3) estimating outcome-specific regression curves. This article proposes a Bayesian approach for addressing this problem, motivated by multisite tumor data from carcinogenicity experiments. A multivariate smoothing spline model is specified, that accommodates dependency in the multiple curves through a hierarchical Markov random field prior for the basis coefficients, while also allowing for residual correlation. A Gibbs sampler is proposed for posterior computation, and the approach is applied to data on body weight and tumor occurrence.