A global/local method with the modified sub-modeling approach of the two-parametric optimal equivalent volume solder balls is introduced to predict the deformation and reliability of the package. The equivalent solder balls as exhibit in this method can obviously reduce the required elements/nodes quantities to enhance computing efficiency. A package model of wire-bonded stacked chip ball grid array under cyclic thermal loading is used as a test vehicle to verify the influences of design factors by fatigue life indicator. Comparing the proposed method with the global fine mesh model, it is found that the difference in the accumulated strain energy density is merely 5.77%, but the optimal equivalent model has highly saved 90% finite element analysis required elements, which means the adopted method can effectively replace the global fine mesh model because both results are in accordance with each other. Using design of experiments to efficiently verify each factor influence with their crosscoupling effects, this paper adopts two kinds of response surface methods that confirm the fatigue life of the proposed approach can be improved by as much as 123.7% for the dual response surface method and 126.3% for the mixed response surface method when comparing with the baseline model. In addition, the optimization of generic algorithm for both response surface methods is demonstrated in this study. From the reviews of factor coupling effects, it is concluded that the response surface method is eligible to achieve the optimum design for package reliability improving.