In this study, we apply an empirical scoring method to evaluate the feasibility of probabilistic seismic hazard analyses at regional scale in Italy accounting for site amplification, which is taken into account through the application of a set of ground motion prediction equations (GMPEs) defined for specific ground types. Precisely, this method calculates the agreement (in terms of likelihood) between the hazard results computed using a specific hazard model and the number of ground-motion exceedances at a set of reference sites. Such a procedure is applied to quantify the likelihood of the outcomes of different hazard models, each based on a specific GMPE, with respect to the observations at 56 accelerometric sites operating in Italy for at least 25 years. Indirectly, this allows evaluating the influence of the selected GMPE in providing reliable hazard estimates. Seven possible GMPEs, applicable in active shallow crustal regions like Italy, have been examined taking into account the correlation among the hazard estimates provided by each computational model at the reference sites. Our results indicate that, although using GMPEs for pre-defined soil categories provides only a broad assessment of the hazard (since it ignores specific site response), large-scale hazard maps of Italy that are compatible with observations can be provided when suitable GMPEs are considered. A restricted number of GMPEs was found appropriate to this scope.