This research seeks to further develop and implement a class of versatile multivariate nonpararnetric statistical models known as the smoothing spline analysis of variance (SSANOVA) models. The research consists of three closely related aspects: algorithms, methodologies, and software implementation. For algorithms, the focus is on the development of a strategy that uses computationally less expensive approximations to the exact smoothing splines: one may reduce the operation from a cubic order to slightly over a linear order while maintaining the same asymptotic convergence rates. For methodologies, empirical studies of some recently proposed (conditional) density and hazard models will be conducted with the assistance of the scalable algorithms, and modeling tools such as diagnostics for the practical significance of model terms will be developed. Among the models to be studied are models for nonparametric regression with arbitrary multiple responses and hazard models more general than the proportional hazard models. Also to be explored is an approach to random effects that allows straightforward numerical treatments for regression (Gaussian and non Gaussian) with correlated errors and for hazard models with frailty. For software implementation, friendly user interface will be implemented in R with the syntax largely inherited from the standard parametric modeling tools. The final product of this research will be a comprehensive R package freely downloadable, with on-line documentations and a separate manual filled with examples. Earlier developments of SSANOVA models in regression have found applications in large epidemiological studies and the modeling of rhythms data and hormones data. The user friendly, scalable software will facilitate the adoption and routine use of these versatile modeling tools by biomedical researchers from a broader spectrum of application areas.