This work presents the random sampling - high dimensional model representation (RS-HDMR) algorithm for identifying complex bionetwork structures from multivariate data. RS-HDMR describes network interactions through a hierarchy of input-output (IO) functions of increasing dimensionality. Sensitivity analysis based on the calculated RS-HDMR component functions provides a statistically interpretable measure of network interaction strength, and can be used to efficiently infer network structure. Advantages of RS-HDMR include the ability to capture nonlinear and cooperative realtionships among network components, the ability to handle both continuous and discrete relationships, the ability to be used as a high-dimensional IO model for quantitative property prediction, and favorable scalability with respect to the number of variables. To demonstrate, RS-HDMR was applied to experimental data measuring the single-cell response of a protein-protein signaling network to various perturbations. The resultant analysis identified the network structure comparable to that reported in the literature and to the results from a previous Bayesian network (BN) analysis. The IO model also revealed several nonlinear feedback and cooperative mechanisms that were unidentified through BN analysis.