Physical contact established during interaction between a human operator and a haptic device creates a coupled system with stability and performance characteristics different than its individual subsystems taken in isolation. Proper incorporation of operator dynamics in physical human-robot interaction (pHRI) conditions requires knowledge of system variables and parameters, some of which are not directly measurable. Operator endpoint impedance, for instance, cannot be directly measured in typical haptic control conditions. Several endpoint impedance estimation techniques have been explored in previous literature, based on measured kinematics and/or other correlated metrics. Arm muscle activity, measured through surface electromyography (sEMG), has been used in previous literature to estimate endpoint stiffness, which is the static component of impedance. Co-activation (co-contraction) of antagonistic arm muscles forming a pair around a joint is known to be the driving factor in modulation of endpoint stiffness. However, previous work employing muscle co-contraction to predict endpoint stiffness has mainly been absent, due in part to the inefficacy of operator models to incorporate muscle co-activation into the prediction scheme. The current study proposes a method for prediction of operator endpoint stiffness based on measured co-contraction levels of a select group of muscles. The proposed methodology incorporates an upper extremity musculoskeletal model that accounts for muscle redundancy and the role of muscle co-contraction on arm stiffness modulation. The study hypothesizes that a free parameter, currently known in the literature to represent the nullspace of the mapping between muscle forces and joint torques, is a random variable of which probability density function can be estimated. Changes in this parameter directly affect changes in operator endpoint stiffness. Ten healthy subjects were asked to resist perturbations induced by a one degree of freedom (DOF) haptic paddle device while measurements, including muscle activities of four arm muscles, were being carried out. Direct derivation of stiffness values at the endpoint was compared to simulated endpoint stiffness values obtained using the proposed predictive methodology. Ten out of forty prediction trials resulted in a statistically significant correlation between predicted and actual stiffness values. Impressive stiffness prediction results, with over 99% peak accuracy in value, were found in only one trial using a combination of the proposed method along with a standard static optimization method for muscle force computation. Though the possibility of using a probabilistic approach to stiffness prediction was shown, robustness and generalizability of the proposed approach to multi-DOF systems remain to be addressed.