The paper focuses on the enhanced stiffness modeling of robotic manipulators by taking into account influence of the external force/torque acting upon the end point. It implements the virtual joint technique that describes the compliance of manipulator elements by a set of localized six-dimensional springs separated by rigid links and perfect joints. In contrast to the conventional formulation, which is valid for the unloaded mode and small displacements, the proposed approach implicitly assumes that the loading leads to the non-negligible changes of the manipulator posture and corresponding amendment of the Jacobian. The developed numerical technique allows computing the static equilibrium and relevant force/torque reaction of the manipulator for any given displacement of the end-effector. This enables designer detecting essentially nonlinear effects in elastic behavior of manipulator, similar to the buckling of beam elements. It is also proposed the linearization procedure that is based on the inversion of the dedicated matrix composed of the stiffness parameters of the virtual springs and the Jacobians/Hessians of the active and passive joints. The developed technique is illustrated by an application example that deals with the stiffness analysis of a parallel manipulator of the Orthoglide family.