  In this paper, we propose and study construction of confidence bands for
shape-constrained regression functions when the predictor is multivariate. In
particular, we consider the continuous multidimensional white noise model given
by $d Y(\mathbf{t}) = n^{1/2} f(\mathbf{t}) \,d\mathbf{t} + d W(\mathbf{t})$,
where $Y$ is the observed stochastic process on $[0,1]^d$ ($d\ge 1$), $W$ is
the standard Brownian sheet on $[0,1]^d$, and $f$ is the unknown function of
interest assumed to belong to a (shape-constrained) function class, e.g.,
coordinate-wise monotone functions or convex functions. The constructed
confidence bands are based on local kernel averaging with bandwidth chosen
automatically via a multivariate multiscale statistic. The confidence bands
have guaranteed coverage for every $n$ and for every member of the underlying
function class. Under monotonicity/convexity constraints on $f$, the proposed
confidence bands automatically adapt (in terms of width) to the global and
local (H\"{o}lder) smoothness and intrinsic dimensionality of the unknown $f$;
the bands are also shown to be optimal in a certain sense. These bands have
(almost) parametric ($n^{-1/2}$) widths when the underlying function has
``low-complexity'' (e.g., piecewise constant/affine).
