High dimensional model representation is under active development as a set of quantitative model assessment and analysis tools for capturing high-dimensional input-output system behavior based on a hierarchy of functions of increasing dimensions. The HDMR component functions are optimally constructed from zeroth order to higher orders step-by-step. This paper extends the definitions of HDMR component functions to systems whose input variables may not be independent. The orthogonality of the higher order terms with respect to the lower order ones guarantees the best improvement in accuracy for the higher order approximations. Therefore, the HDMR component functions are constructed to be mutually orthogonal. The RS-HDMR component functions are efficiently constructed from randomly sampled input-output data. The previous introduction of polynomial approximations for the component functions violates the strictly desirable orthogonality properties. In this paper, new orthonormal polynomial approximation formulas for the RS-HDMR component functions are presented that preserve the orthogonality property. An integrated exposure and dose model as well as ionospheric electron density determined from measured ionosonde data are used as test cases, which show that the new method has better accuracy than the prior one.