This contribution investigates potential best practice in statistical treatment of exploration geochemistry data for drill target definition, using data from the Sari Gunay epithermal gold deposit, northwest Iran. Primarily, a robust factor analysis was applied to the soil geochemical data to identify anomalies related to outcropping and weak surface expression of hidden gold mineralizations. Investigation of outcropping mineralization indicates that the invisible gold mostly occurs in arsenian pyrite, accompanied by antimony and mercury sulfides. Therefore, a favorable principal factor with strong association of Au, As, Sb, Hg, Tl and S elements, related to epithermal mineralization, was selected for analysis. Then, the boxplot and median. +. 2MAD (median absolute deviation) techniques of the exploratory data analysis (EDA) and also concentration-area (C-A) fractal model were applied to the favorable factor loading to set the threshold values for defining multi-element soil geochemical anomalies for drilling. The median. +. 2MAD and C-A fractal models successfully identified the known outcropping gold mineralization at the Sari Gunay and Agh Dagh hills. In contrast, the boxplot method only identified the strong outcropping gold mineralization at the Sari Gunay hill, but failed to map the weaker mineralization at the Agh Dagh hill due to the usage of higher threshold values. In addition, the C-A fractal analysis identified a few second-class, weaker multi-element soil anomalies, mostly related to subsurface gold mineralization. The identified soil geochemical anomalies were compared with the in situ economic gold concentration in channel rock samples, mostly below the soil cover, and the outcropping host lithologies for validation. The results indicated the EDA approaches failed to map weak surface expression of gold mineralization and only mapped strong geochemical anomalies, useful only for reconnaissance drilling stages. The C-A fractal analysis is capable of mapping multi-class geochemical anomalies, some of which were associated with the weak surface expression of hidden gold mineralization, and is therefore useful in defining targets at both the reconnaissance and systematic drilling stages. The identified geochemical anomalies by C-A fractal analysis were also useful to characterize the favorable lithologies (e.g. quartz-tourmaline breccias, intrusive fragmental rocks including diatreme and dacite/trachyte porphyry) and hydrothermal alterations (e.g. quartz-sericite, K-feldspar and silicification) hosting gold mineralization at Sari Gunay. © 2014 Elsevier B.V.