Eighteen duplicate, composite soil samples from heterogeneous remote historic metal ore mining sites (Miedzianka Mt. and Karczowka Mt., south-central Poland) were analyzed twice for As, Cd, Co, Cr, Cu, Mn, Ni, Pb and Zn by the ICP-MS method after aqua regia extraction. Subsequently, the results were tested for normality. The sampling uncertainty [expressed as the relative standard deviation srsamp (%)] was computed using three different methods: ANOVA, RANOVA and range statistics. However, it was possible to use all three methods only for Cr and Pb (Miedzianka) and for Co (Karczowka). The obtained values were in the ranges of: 15.6–16.0 % for Cr, 17.0–20.2 % for Co and 24.6–38.7 % for Pb. The sampling uncertainty for these elements that revealed non-normally distributed data was calculated using the robust ANOVA. The results varied from 12.1 % (As) to 31.9 % (Ni) in soils from Karczowka, and from 9.2 % (Co) to 35.9 % (Cu) in soils from Miedzianka. The highest srsamp values exceeding 20 % were noted for Cd (31.9 %) and Ni (26.0 %) in soils from Karczowka, and for Cu (20.5 %) and Mn (35.9 %) in soils from Miedzianka. Because traditional methods of transformation were insufficient to reduce non-normality, the Box–Cox model was applied to these four elements. The sampling uncertainty, computed for transformed data with a one-way ANOVA method, was as follows: 20.1 % (Cd), 7.3 % (Ni), 4.6 % (Cu) and 4.6 % (Mn) and for back-transformed data: 39.7 % (Cd), 12.8 % (Ni), 7.4 % (Cu) and 12.4 % (Mn). The sampling uncertainty was lower compared to the values calculated with raw data, but the interpretation of results obtained was problematic.