System identification as a dynamic system modelling approach can, in principle, be applied to determine the dynamic model of building structures during earthquakes and conduct rapid building damage assessment via time or frequency analysis of the identified building model. However, under large-amplitude earthquake ground motions, building structures behave nonlinearly, making well-established linear system identification and associated analysis not applicable in many cases. To address this issue, many studies have been conducted to take structural nonlinearity into account while performing system identification-based building damage assessment. However, so far, only specific nonlinearity such as stiffness nonlinearity has been considered. There are still no approaches that can be applied to a wide range of nonlinear building damage assessment problems. In order to resolve this challenge, in the present study, a data driven nonlinear modelling and model frequency analysis approach is proposed. The approach applies a nonlinear system identification technique to determine the Nonlinear Autoregressive with eXogenous input (NARX) model of building structures during earthquakes using recorded ground motion data and corresponding building acceleration responses. Then the nonlinear output frequency response functions (NOFRFs) of the NARX model are used to assess the building damage induced by earthquakes. Simulation studies are conducted where the earthquake responses of an elastic-plastic mass-spring-damping model of building structures are used to verify the effectiveness of the proposed NARX modelling and model NOFRFs analysis approach. After that, the proposed approach is applied to the E-Defense shaking table data to demonstrate its practical significance and potential engineering applications.