A Robust algorithm for the extraction of reduced-order behavioral models from sampled frequency responses is proposed. The system under investigation can be any Linear and Time Invariant structure, although the main emphasis is on devices that are relevant for Signal and Power Integrity and RF design, such as electrical interconnects and integrated passive components. We assume that the device under modeling is parameterized by one or more design variables, which can be related to geometry or materials. Therefore, we seek for multivariate macromodels that reproduce the dynamic behavior over a predefined frequency band, with an explicit embedded dependence of the model equations on these external parameters. Such parameterized macromodels may be used to construct component libraries and prove very useful in fast system-level numerical simulations in time or frequency domain, including optimization, what-if, and sensitivity analysis. The main novel contribution is the formulation of a finite set of convex constraints that are applied during model identification, which provide sufficient conditions for uniform model stability and passivity throughout the parameter space. Such constraints are characterized by an explicit control allowing for a trade-off between model accuracy and runtime, thanks to some special properties of Bernstein polynomials. In summary, we solve the longstanding problem of multivariate stability and passivity enforcement in data-driven model order reduction, which insofar has been tackled only via either overconservative or heuristic and possibly unreliable methods.