The restoring force method (RFM), a nonparametric identification technique established in applied mechanics, was used to maximize the information obtained from moment-rotation hysteresis curves under pure moment flexion-extension testing of human lumbar spines. Data from a previous study in which functional spine units were tested intact, following simulated disk injury, and following implantation with an interspinous process spacer device were used. The RFM was used to estimate a surface map to characterize the dependence of the flexion-extension rotation on applied moment and the resulting axial displacement. This described each spine response as a compact, reduced-order model of the complex underlying nonlinear biomechanical characteristics of the tested specimens. The RFM was applied to two datasets, and successfully estimated the flexion-extension rotation, with error ranging from 3 to 23%. First, one specimen, tested in the intact, injured, and implanted conditions, was analyzed to assess the differences between the three specimen conditions. Second, intact specimens (N = 12) were analyzed to determine the specimen variability under equivalent testing conditions. Due to the complexity and nonlinearity of the hysteretic responses, the mathematical fit of each surface was defined in terms of 16 coefficients, or a bicubic fit, to minimize the identified (estimated) surface fit error. The results of the first analysis indicated large differences in the coefficients for each of the three testing conditions. For example, the coefficient corresponding to the linear stiffness (a01) had varied magnitude among the three conditions. In the second analysis of the 12 intact specimens, there was a large variability in the 12 unique sets of coefficients. Four coefficients, including two interaction terms comprised of both axial displacement and moment, were different from zero (p < 0.05), and provided necessary quantitative information to describe the hysteresis in three dimensions. The results suggest that further work in this area has the potential to supplement typical biomechanical parameters, such as range of motion, stiffness, and neutral zone, and provide a useful tool in diagnostic applications for the reliable detection and quantification of abnormal conditions of the spine.