In this work, the viscoelastic behavior of high-density polyethylene (HDPE) and polypropylene (PP) was studied through stress relaxation experiments conducted at different strain levels. The main objective was to evaluate classical, fractional, and conformable derivatives to analyze molecular mobility, using statistical methods to identify the most accurate representation of the viscoelastic response. Besides the coefficient of determination (R²), the average absolute deviation (AAD) and mean squared error (MSE) were used as evaluation metrics, along with a multivariate analysis of variance (MANOVA) and the response surface methodology (RSM) to optimize the correspondence between experimental data and model predictions. The findings demonstrate that the spring-pot, Fractional Maxwell (FMM), Fractional Voigt–Kelvin (FVKM), and Kohlrausch–Williams-Watts (KWW) models effectively describe stress relaxation under statistical criteria. However, a joint analysis using RSM revealed that the choice of mathematical model significantly influences the outcomes. The FVKM was identified as the most effective for HDPE, while the KWW model best characterized PP. These results highlight the importance of optimization tools in advancing the characterization of polymer viscoelasticity. The ability to select the most accurate models for HDPE and PP under varying conditions can directly improve the performance and durability of products in critical industrial sectors such as packaging, automotive, and medical devices, where long-term mechanical behavior is crucial. By offering a framework adaptable to other materials and modeling approaches, this work provides valuable insights for optimizing polymer processing, improving product design, and enhancing the reliability of polymer-based components in a range of industrial applications.