A nonlinear response surface method incorporating multivariate spline interpolation (RSM-S) is a useful technique for the optimization of pharmaceutical formulations, although the direct reliability of the optimal formulation must be evaluated. In this study, we demonstrated the feasibility of using the bootstrap (BS) resampling technique to evaluate the direct reliability of the optimal liposome formulation predicted by RSM-S. The formulation characteristics ( X _n ), including vesicle size ( X ₁), size distribution ( X ₂), zeta potential ( X ₃), elasticity ( X ₄), drug content ( X ₅), entrapment efficiency ( X ₆), release rate ( X ₇), and the penetration enhancer (PE) factors as formulation factors ( Z _n ), with the type of PE ( Z ₁) and content of PE ( Z ₂) were used as causal factors of the response surface analysis. The intended responses were high skin permeability (flux, Y ₁) and high stability formulation (drug remaining, Y ₂). Based on the dataset obtained, the simultaneous optimal solutions were estimated using RSM-S. Leave-one-out-cross-validation showed satisfying reliability of the optimal solution. Concurrently, similar BS optimal solutions were estimated from the BS dataset that was generated from the original dataset through BS resampling at frequencies of 250, 500, 750, and 1000. The analysis and simulation indicated that X ₄, X ₅, and Z ₂ were the prime factors affecting Y ₁ and Y ₂. These findings suggest that this approach could also be useful for evaluating the reliability of an optimal liposome formulation predicted by RSM-S and would be beneficial for the pharmaceutical development of liposomes for transdermal drug delivery.