Potential earthquake sources are revealed by the locations of earthquakes in historical and instrumental catalogs and of geologically mapped faults, including plate boundaries. We derive a set of multiplicative hybrid earthquake likelihood models that combine earthquake and fault data for the New Zealand region. In these models, the cell rates in a spatially uniform baseline model are scaled using selected subsets of five covariates derived from the magnitudes and locations of past earthquakes, the location of the boundary between the Australian and Pacific plates, and the location and slip rate of mapped faults. The hybrid model parameters are optimized for earthquakes of M 5 and greater over the period 1987–2006 and tested on earthquakes from the period 2007–2014. No updating of models is undertaken during the fitting or testing period, but we consider two cases of the earthquake-based covariates in the tests: (1) all data prior to 1987 and (2) all data prior to 2007, respectively. Hybrids containing the earthquake-based covariates perform better in the latter case. The most informative hybrid models in the fitting and testing period are composed of three and four covariates, respectively, including both earthquake- and fault-based variables. Proximity to mapped faults is overall the most informative individual covariate. These results can be used to inform better modeling of long-term earthquake occurrence rates for probabilistic seismic-hazard analysis.