The present study was conducted to determine different geochemical anomalies of rare earth elements (REEs) using a combined approach of stepwise factor analysis (SFA), sequential Gaussian simulation (SGS) ,and concentration-area (C-A) fractal modeling based on surface lithogeochemical samples obtained from the Esfordi phosphate mine (Central Iran). The Esfordi mine is one of the important mines in Bafq metallogenic zone due to average and maximum grade of 0.5 and 1.7%, respectively for REEs. According to SFA operation in two steps, the REEs were placed in the first factor of the second stage (F1-2). Then, the SGS method and C-A fractal modeling were performed on F1-2 factor scores for classification of anomalies. Log-ratio matrix was used to evaluate the correlation of these results with anomalous lithogeochemical samples, as well as to determine the relationship of anomalies with rock types and mineralized units and finally, to validate the results of SGS-fractal modeling. The results confirmed an appropriate correlation between F1-2 anomalies and high-concentration in further rock samples. The main anomalies were found to have good correlation with apatite-iron unit and in general with other apatite-bearing units based on overall accuracy values. The apatite-bearing units with high values of REEs were located in northern and central parts of the mine. The results of combined approach of SFA, SGS, and C-A fractal modeling showed that this hybrid approach can be useful in determining anomalies with proper accuracy.