This paper presents a systematic and efficient design approach for the two degree-of-freedom (2-DoF) capacitive microelectromechanical systems (MEMS) accelerometer by using combined design and analysis of computer experiments (DACE) and Gaussian process (GP) modelling. Multiple output responses of the MEMS accelerometer including natural frequency, proof mass displacement, pull-in voltage, capacitance change, and Brownian noise equivalent acceleration (BNEA) are optimized simultaneously with respect to the geometric design parameters, environmental conditions, and microfabrication process constraints. The sampling design space is created using DACE based Latin hypercube sampling (LHS) technique and corresponding output responses are obtained using multiphysics coupled field electro–thermal–structural interaction based finite element method (FEM) simulations. The metamodels for the individual output responses are obtained using statistical GP analysis. The developed metamodels not only allowed to analyze the effect of individual design parameters on an output response, but to also study the interaction of the design parameters. An objective function, considering the performance requirements of the MEMS accelerometer, is defined and simultaneous multi-objective optimization of the output responses, with respect to the design parameters, is carried out by using a combined gradient descent algorithm and desirability function approach. The accuracy of the optimization prediction is validated using FEM simulations. The behavioral model of the final optimized MEMS accelerometer design is integrated with the readout electronics in the simulation environment and voltage sensitivity is obtained. The results show that the combined DACE and GP based design methodology can be an efficient technique for the design space exploration and optimization of multiphysics MEMS devices at the design phase of their development cycle.