In this paper, the long-term extreme response of a vessel rolling in random beam seas is addressed. The long-term response analysis is based on the upcrossing rates of the roll motion under different sea states. However, the nonlinear effects associated with the restoring and damping terms have a significant influence on the high-level response, assessing the corresponding statistics, such as the upcrossing rate, with low probability levels is difficult and time-consuming. In this work, the Markov theory is introduced in order to tackle this problem. Specifically, the random roll excitation moment is approximated as a filtered white noise process by applying a linear filter technique and an efficient four-dimensional (4D) path integration (PI) procedure is applied in order to calculate the response statistics. The long-term analysis of nonlinear roll motion in random seas that takes into considerations of the response statistics obtained by the 4D PI method as well as the variation of the sea states could be a valuable reference for ship stability research.