The paper is devoted to the stiffness modeling of heavy industrial robots with pneumatic gravity compensators. The main attention is paid to the identification of elastostatic parameters and calibration accuracy. To identify the desired set of parameters two-step identification procedure is developed. To reduce impact of the measurement errors, the set of manipulator configurations for calibration experiments is optimized with respect to the proposed performance measure related to the end-effector position accuracy. Advantages of the developed technique are illustrated by dedicated experimental study.