  We consider the nonparametric regression problem with multiple predictors and
an additive error, where the regression function is assumed to be
coordinatewise nondecreasing. We propose a Bayesian approach to make an
inference on the multivariate monotone regression function, obtain the
posterior contraction rate, and construct a universally consistent Bayesian
testing procedure for multivariate monotonicity. To facilitate posterior
analysis, we set aside the shape restrictions temporarily, and endow a prior on
blockwise constant regression functions with heights independently normally
distributed. The unknown variance of the error term is either estimated by the
marginal maximum likelihood estimate or is equipped with an inverse-gamma
prior. Then the unrestricted block heights are a posteriori also independently
normally distributed given the error variance, by conjugacy. To comply with the
shape restrictions, we project samples from the unrestricted posterior onto the
class of multivariate monotone functions, inducing the "projection-posterior
distribution", to be used for making an inference. Under an
$\mathbb{L}_1$-metric, we show that the projection-posterior based on $n$
independent samples contracts around the true monotone regression function at
the optimal rate $n^{-1/(2+d)}$. Then we construct a Bayesian test for
multivariate monotonicity based on the posterior probability of a shrinking
neighborhood of the class of multivariate monotone functions. We show that the
test is universally consistent, that is, the level of the Bayesian test goes to
zero, and the power at any fixed alternative goes to one. Moreover, we show
that for a smooth alternative function, power goes to one as long as its
distance to the class of multivariate monotone functions is at least of the
order of the estimation error for a smooth function.
