The Cox point process is highly considered for earthquake modeling. However, the complex earthquake data which involve a large number of occurrences and geological variables often require expensive computation. This study aims to propose an efficient algorithm based on the two-step procedure by constructing the first and second order composite likelihoods. We consider four Neyman–Scott Cox process models and apply them to fit the earthquake distribution in Sumatra. We conclude that the Cauchy cluster process performs best.