One of important factors in designing MEMS resonators for RF filters is obtaining a desired frequency response function (FRF) within a specific frequency range of interest. Because various array-type MEMS resonators have been recently introduced to improve the filter characteristics such as bandwidth, pass-band, and shape factor, the degrees of freedom (DOF) of finite elements for their FRF calculation dramatically increases and therefore raises computational difficulties. In this paper the Krylov subspace-based model order reduction using moment-matching with non-zero expansion points is represented as a numerical solution to perform the frequency response analyses of those array-type MEMS resonators in an efficient way. By matching moments at a frequency around the specific operation range of the array-type resonators, the required FRF can be efficiently calculated regardless of their operating frequency from significantly reduced systems. In addition, because of the characteristics of the moment-matching method, a minimal order of reduced system with a prearranged accuracy can be determined through an error indicator using successive reduced models, which is very useful to automate the order reduction process and FRF calculation for structural optimization iterations. We also found out that the presented method could obtain the FRF of a $6 imes6$ array-type resonator within a seventieth of the computational time necessary for the direct method and in addition FRF calculation by the mode superposition method could not even be completed because of a data overflow with a half after calculation of 9,722 eigenmodes.