Site response models are frequently used in engineering practice to predict surficial ground motions based on a site-specific soil profile and input motions, and site response predictions are especially important for large strains and accelerations, which have a greater damage potential. To characterize nonlinear soil behavior at large strains, a number of constitutive soil models have been developed. However, the application of fully nonlinear time-domain site response analyses remains limited in practice, with the equivalent-linear site response approximation to nonlinear soil behavior, using frequency-domain programs such as SHAKE (Schnabel et al., 1972), still the most common approach. For a particular project, engineering practitioners are therefore faced with the challenge of selecting the appropriate level of model complexity (e.g., equivalent-linear vs. nonlinear). While previous validation studies have attempted to quantify the levels of ground motion for which nonlinear site response analyses are necessary (e.g., Assimaki et al., 2008; Kwok et al., 2008; Kim and Hashash, 2013; Kaklamanos et al., 2015), the assessment of fully nonlinear site response models is often limited to a relatively small number of sites and ground motions. In this study, one-dimensional (1D) total-stress nonlinear, equivalent linear, and linear site response predictions are calculated using an unprecedented number of sites and ground motions, allowing for more statistically significant conclusions to be drawn than in prior studies. This study uses Japan’s comprehensive Kiban-Kyoshin network of vertical seismometer arrays (Aoi et al., 2000), in particular, 5626 ground motions at 114 KiK-net sites are utilized, with 239 ground motions having PGA > 0.3g. Site response predictions are calculated using the program DEEPSOIL (Hashash et al., 2011), and SHAKE for the nonlinear, and equivalent linear analyses, respectively; based on the P- and S- wave velocity profiles, and soil types provided on the KiK-Net database. The Zhang et al. (2005) modulus-reduction and damping curves are used in the equivalent-linear analyses and as the target curves for the nonlinear analyses. This study builds upon prior work (Kaklamanos et al., 2013) in which linear and equivalentlinear site response analyses (but not nonlinear analyses) were performed at 100 KiK-net sites using 3720 ground motions, allowing for broad conclusions on the uncertainty of linear and equivalent-linear site response models. With the large database of nonlinear site response model predictions in the current study, the predictive capabilities of fully nonlinear total-stress site response models relative to linear and equivalent-linear models are assessed. The model residuals assessed in this study are those of the 5%-damped pseudo-acceleration response spectra, calculated as ln(PSAobs) – ln(PSApred), where PSAobs and PSApred are the observed and predicted spectral accelerations at a given period, respectively. From analyzing the trends of the model residuals versus the maximum shear strain in the soil profile, Kaklamanos et al. (2013) concluded that the equivalent-linear model becomes inaccurate when strains exceed 0.1 to 0.4%. In the current study, we find that the model residuals of the equivalent-linear and nonlinear site response models generally do not deviate from each other significantly at large shear strains. For shear strains greater than 0.5% at short spectral periods, both the equivalent-linear and nonlinear model residual plots slope upwards, indicating that these models tend to underpredict large-strain ground motions. However, the nonlinear model residuals do not slope upward as significantly at some spectral periods (for example, for spectral 1 accelerations at T = 0.1 s). Furthermore, the scatter in the equivalent-linear model residuals is greater than that of the nonlinear model residuals at large shear strains, suggesting that the equivalent-linear site response model is less precise at large shear strains. In the aggregate, the linear, equivalent-linear, and nonlinear model biases and standard deviations can be calculated across all sites and ground motions using mixed-effects regression on the model residuals. Comparisons of the model biases and standard deviations indicate that all 1D site response models (linear, equivalent-linear, and nonlinear) are biased towards underprediction of ground motions at short spectral periods, where nonlinear effects are strongest. However, the equivalent-linear and nonlinear model biases are smaller than the linear model bias. The persistent model biases suggest that: (1) many of these sites may experience a breakdown in the 1D site-response assumptions; and/or (2) the site investigation data provided on KiK-net (i.e. velocity profiles and broad soil type) may be over-simplified. With respect to the first point, in particular, the underlying assumptions of 1D site response may have to be addressed in order to make notable prediction improvements, perhaps by incorporation of three-dimensional soil constitutive response and incident ground motion effects. Based on the inter-site residuals, we have also identified some “interesting” sites at which all 1D site response models most strongly overpredict or underpredict ground motions: ISKH05 and KOCH05 are characterized by the strongest underpredictions, and HYGH07, IWTH07, and WKYH01 are characterized by the strongest overpredictions (at different vibration periods, however). Because these site-specific biases are consistent across all 1D site response models, the 1D site response assumption is likely not valid at these sites. Although the nonlinear site response models are shown to offer an improvement over equivalent-linear models, the remaining trends in the nonlinear model residuals suggest that other factors—such as three-dimensional effects—have a significant impact on site response behavior.