Thursday, 27 April 2017

Extraordinary claims require extraordinary evidence

The history of claims for the earliest human presence in the Americas is littered with extraordinary claims, many of which have not stood up to scrutiny. In many cases it has been the dating rather than the human associations of the site that has failed to withstand other scientists' examination. This week Nature has published a paper by Holen and coworkers on “A 130,000-year-old archaeological site in southern California, USA”. Of great interest to me is that the dating is primarily based on “230Th/U radiometric analysis of multiple bone specimens using diffusion–adsorption–decay dating models”. The first development of this approach to dating bone was in my doctoral thesis and it was published in my first scientific paper. But I've also published quite a few papers on re-evaluating the claimed dates for Palaeolithic sites and the minimum level of information that needs to be published for readers to be able to assure themselves of the quality of dates. So, how does this look?

The uranium-series dates

The authors are to be congratulated in making the full set of radiometric data available in the supplementary information, and making it available in machine-readable from from the USGS ScienceBase. The main text reports a series of checks on the quality of the analytical chemistry, and these are all passed with flying colours. The key to understanding the quality of the uranium-series dates is in Extended Data Figure 10 (this is free to access even though the main paper is not). The diffusion-adsorption model makes a series of assumptions and predictions that are not directly involved in the date calculation using the updated approach used here. The things to watch out for are:
  1. the concentrations of uranium through the bone should be U-shaped, high at the edges and lower in the middle. In the top panel of part (a) of the figure this appears to be broadly the case for all three bones dated, but CM-292 has a decrease towards the right of the plot, which might indicate inhomogeneity in the bone or leaching of uranium (which the model assumes does not happen).
  2. the distribution of apparent uranium-series ages calculated without making any correction for uranium uptake, as shown in the lower panel of part (a) of the figure should also be U-shaped. For this test the data is noisier but follows the expected trends. Here CM-292 is the best of the three samples. The iDAD software used to calculated ages can accommodate scatter like this – it makes the final date less certain.
  3. comparison of the predicted isotope ratios from the best fit model to the actual data, as shown in Panel (c) in the second and third columns. CM-292 clearly has the best fit, and fairly smoothly varying observations, while the other two samples have more jagged distributions of the observations and poorer fits. This is then reflected in,
  4. comparison of dates and their uncertainties calculated via two statistical methods using the same model, shown in the first column of panel (c). If the two approaches agree then we can have confidence in the robustness of the results. This is not well explained in the paper and I had to go back to the paper describing the diffusion-adsorption-decay model to understand what is shown The smooth blue curves are calculated using a Bayesian approach which evaluates the whole probability distribution, while the red lines use a maximum-likelihood approach and assume a normal distribution. The paper does not give ranges based on the Bayesian approach, but it looks to me that for CM-292 the two results are almost the same. For the other two bones the two approaches differ, both in the most likely date, which is younger in the Bayesian approach, and in the uncertainty range, which appears to be wider in the Bayesian approach. But all-in-all the values are fairly close, and  the differences in best-estimates are of no consequence given the uncertainties.
  5. From the data it is possible to calculate the ratio of 234U/238U that was supplied to the bone from groundwater before any radioactive decay took place in the bone. Comparing this to expected values can help test whether the model is correct. The authors state “Initial 234U/238U activity ratios calculated for bone subsamples span a narrow range (1.38-–1.50) that is consistent with modern shallow groundwater from the nearby Sweetwater River drainage (1.45–1.54)”. They do not provide any argument to show that these waters should be comparable to those the bones have been exposed to over a full glacial cycle, nor any indication of how 'nearby' this is. Nevertheless, if this is the right comparison, the closeness of the values is encouraging.
Overall the authors have carried out almost all the checks that could be done using uranium-series data, though they have not evaluated other elements in the bone which could have added to the robustness of the results. In terms of uranium-series dating of bone, this is one of the best-evidenced studies to date.

Luminescence dating

The supplementary material reports a series of luminescence dating studies by Steven Forman at University of Illinois at Chicago and Ronald J. Goble at the University of Nebraska. Neither of these contributors are authors of the article. These attempts were unsuccessful as the upper limit of dating at this site is about 70,000 years. It would have been nice to see more detail of why this upper limit is reached, including plots showing how close to saturation the dates are.

Stratigraphic context

Scientific dating methods don't exist in a vacuum. The dated samples come from a stratigraphic context, which can provide constraints and checks on the dates. The supplementary information describes the site as lying “within a 12 m-thick series of flat-lying, late Pleistocene fluvial sediments”, but the figures and descriptions of sediments only cover about 3.5m of that sequence. To check the dates it would be useful to have a figure of the wider stratigraphic context, showing stratigraphic correlations with nearby sites and previous estimates of the ages of the strata. Without this it is difficult to put what is effectively a single high-quality date into context.

The supplementary material does refer to a site report which may contain more detail. I've not been able to find this online but I did find a scan of the executive summary. This does not contain any further stratigraphic information but does state “Radiometric dating of ivory and soil carbonate from the quarry yielded dates of 335+/-35 Ka (thousands of years before present) and 196+/-15 Ka”. What are these dates? They do not feature in the paper or the supplementary information, but prima facie they differ greatly from the date proposed in the paper and it would seem that this radiometric data needs to be reported and shown to be wrong before we can be confident about the date of the site.


The Cerutti Mastodon site has one well-dated mastodon, but the lack of stratigraphic context and the hints of unpublished contradictory evidence weaken the robustness of the claimed date. I'm not so well placed to evaluate the evidence for human activity at the site, but the lack of formal stone tools, lack of butchery marks and selection of bones and teeth to smash do not look to be like other sites of this period. Further confirmation will be needed before this enters textbooks as the earliest human occupation in the Americas.

Friday, 25 November 2016

A look at the GPS of Elhaik et al.

The GPS (Geographic Population Structure) method is claimed to be a method that can predict the location of an individual's ancestors 1000 years ago, at least according to the University press release that announced it. I'm interested in ancient migrations and statistics in archaeology, and I follow the application of genetics in archaeology. So it sounds like just my thing!

The GPS was first published by Elhaik et al. (2014) and, despite the authors’ disclaimer of financial interests, immediately formed the basis of a commercial offering at Prosapia Genetics. Recently competing financial interests were acknowledged in a corrigendum to the original paper published two and a half years after the original despite the conflict being obvious within a few days of publication. Since publication, the method has been the subject of some online criticism and has recently been offered by another company, GPS Origins. Critical overviews of these offerings have been given by Debbie Kennett in her blog and Pavel Flegontov on Facebook. GPS has formed the basis for some subsequent publications from the same group. One on the origins of the Jews (Das et al. 2016) has been notably controversial (Flegontov et al.2016, Aptroot  et al. 2016, and see the set of links provided by Pavel Flegontov on Facebook). Recently Das et al. published a preprint of a lengthy response to their critics (Das et al 2016).

Prompted by this I have gone back to look at the original paper to understand its methods a bit better.

The method is simple. Elhaik et al. take genetic data from a set of individuals of known population origin (I shall call this the reference dataset), and compute ADMIXTURE proportions of putative ancestral populations for those individuals. Comparison of the differences in the proportions to the geographical distance between the populations reveals a linear relationship. This relationship is then exploited to locate genetically characterised individuals of unknown geographical origin (I shall call this the test sample) with respect to the reference populations.

However, the devil is in the detail. Joe Pickrell has blogged about the conceptual problems of relating genetic variation to spatial distance in this way. I’m not going to comment on the genetics, but I do know about prediction of spatial locations from other data. The rest of this post covers some mathematical and technical aspects of the paper. For those who wish to cut to the chase, I give a summary and conclusions here.

The regression model treats composition data as if it did not have to sum to 100% which is clearly untrue. The results in the paper cannot be replicated from the supplementary material, using either the method described in the paper or the one actually implemented in the code. The paper describes removing points from the dataset but unless the removed points are improperly replaced with zeroes, the correlation is too weak to be useful for prediction. Consequently, the regression equation calculated is probably invalid as it is methodologically unsound. The supplementary code does not do what is described in the paper.

The predictive model implemented in the code does not use the equations in the paper. In addition, in both the paper and the code if there were two populations of the same genetic distance from the test sample, but different geographical distance from the population that best matches the test sample, the geographically more distant population will contribute more to determining the direction of the GPS arrow. Usually more distant populations should contribute less to a prediction not more.

Test calculations: the (previously available) online tool and downloadable code differ by hundreds of miles in their predicted location for the same data.

Overall conclusion: the mathematical methods described are incoherent, the supplementary data is not that used to create the figures or equations in the paper, and the supplementary code does not implement the methods described. The paper is methodologically unsound and not reproducible.

Now here's the detail ...

The regression model

Elhaik et al. (2014) calculated the genetic distances between individuals as “Euclidean distances between … admixture proportions for the analysed individuals”. This is a mistake. Euclidean distances, are designed to apply to unconstrained variables which can take any value between minus infinity and plus infinity. With compositional data, including admixture proportions, the values are not unconstrained as they must sum to 1 or 100%. It has been known since at least 1897 that applying Euclidean methods to compositional data leads to errors, particularly overestimation of correlation coefficients (Pearson 1897, Aitchison 1986, Aitchison 2005).

Elhaik et al. (2014) give a confusing account of how they calculated geographical distances. They state “we calculated the Euclidean distances between all the populations within the GEN and GEO data sets” but then “distances ΔGEO(X, Y) between two individuals with known latitude and longitude were calculated using the Haversine formula”. The Haversine formula is not Euclidean; it is a method for calculating great circle distances on the earth’s surface (as explained on Wikipedia). This is much better than a Euclidean distance. However, the publically released R program code reads latitude and longitude data and then applies the default method of R's dist function. That default method is Euclidean distance! Given that the distance represented by a degree of longitude varies with latitude, this measure is not one that is even approximately linearly related to distance in miles. It is impossible that Figure 9 of Elhaik et al. could have been produced using the code that has been released to accompany the paper.

The third step is to correlate the (incoherently calculated) genetic distances and geographical distances. Each individual in the reference dataset is compared to all the others and the points are analysed. But there is a big problem here. This approach overstates the number of independent data points. If I know that A is close to both B and C, then I know a lot about the possible values of distance between B and C before I measure it. The distance B-C is not completely independent of the distances A-B and A-C. Regression methods assume independence of all the data points, and if the points are not independent the reliability of the regression is over-estimated.

The paper notes that the relationship is linear only up to about 5000 miles geographical distance so they remove datapoints above 4000 miles distance. The published code does something very strange at this point. Having calculated the distances, it replaces any point with geographical distance over 70 or genetic distance over 0.8 with a point (0,0) and then runs the regression. 70 degrees of latitude is approximately 4800 miles, so this may be an approximation to the 4000 mile cut-off described in the paper. But removing genetic distances over 0.8 is not mentioned in the paper. These points are not removed from the regression in the manner described in the paper; they are moved to lie at the origin. As these are 67% of the data points in the supplementary data files (on as of 16 August 2016), this very strongly weights the regression to pass close to the origin. As coded, the regression line obtained is ΔGEO = 0.94(±0.09) + 68.2(±0.4) ΔGEN, with r2 = 0.78; leaving out these points it is ΔGEO = 9.50(±0.49) + 51.8(±1.1) ΔGEN with r2 = 0.41. This latter equation is rather awkward: it would predict that even at zero genetic distance an average geographical distance of 9.5 “degrees” is expected.

If the data provided with the code is converted to miles using the Haversine formula it is clearly a different dataset to that in the paper (see figures below). As the first is stated to be of individuals and the second of reference populations this is not very surprising.

figure 9 from Elhaik et al., described as plotting individuals, 
reference population data provided with the code, converted to miles using the Haversine formula.

The paper reports a regression equation of ΔGEO = 38.7 + 2523 ΔGEN. Running the regression on the supplementary dataset converted to miles but limited to geographical distances less than 4000 miles gives: ΔGEO = 747(±22) + 2422(±35) ΔGEN (r2 = 0.50) which is not very similar to the equation reported in the paper. Removing the points with ΔGEN>0.8 gives ΔGEO = 504(±27) + 3095(±58) ΔGEN (r2 = 0.46). Replacing the removed points (67% of the data) with (0,0) gives ΔGEO = 50(±5) + 3969(±19) ΔGEN (r2 = 0.81).

Summary: the results cannot be replicated from the supplementary material, using either the method described in the paper or the one actually implemented in the code. Unless the removed points are improperly replaced with (0,0), the r2 value is unacceptably low for prediction (below 50%).

Conclusion: the regression equation calculated is probably invalid as it is methodologically unsound The supplementary code does not do what is described in the paper.

The predictive model

Having established their regression equation, Elhaik et al. propose a method to predict the location of a new genetic sample. They calculate the genetic distance to the closest population and use the regression equation to calibrate it to a geographical distance. The sample is assumed to lie at that distance from the closest population and the next nine nearest populations are used to determine the direction of that distance by calculating a weighted average of their vectors from the best match. How the admixture proportions attributed to populations are calculated from the data for individuals is not stated. Presumably they are means. Why a regression equation calculated from distances between pairs of individuals can be applied to the distance between an individual and a population mean is not explained.

In the code the ΔGEOmin value is calculated not with the full regression equation but only by multiplying the ΔGEOmin by the gradient; the intercept is omitted!

The weights for the second to tenth distances given in the paper are

but in the code these are given by

W <- (minE[1]/minE)^4;
W =  W/(sum(W));

As the minE of the code are the ΔGEN of the equations in the paper, the ratio in the first line makes sense, but somehow a fourth power has appeared and then a normalisation so that weights sum to one. Bizarrely the way the code and the paper are constructed, if there were two populations of the same genetic distance but different geographical distance, the geographically more distant population will contribute more to determining the direction of the result because its vector will be longer.

The equations in the paper are constructed so that the result is a distance ΔGEOmin from the best matching population. The code does not do this. The code places the result at that distance or closer if the weighted average vector is shorter.

Conclusion: the code does not implement the predictive equations of the paper, and in any case the method using inappropriate weightings..

Test calculations

The download site also provides an online calculator and partial results for the test data. Entering the first line of test data in the online calculator gives a predictions of Latitude: 34.4237, Longitude: 52.7614 (as of 16 August 2016: The partial results file has 39.5307070554391, 46.885015601252 (yes the values are reported to a precision better than the size of an atom), and is the same as running the supplied code. The instructions on the download page ( state that the online tool gives Latitude: 38.2188 Longitude: 47.2863 for that first line.  [Note: as of 25 November 2016 this site seems to have been hacked and now redirects to a Chinese site.]

Conclusion: the online tool and downloadable code differ by hundreds of miles in the predicted location.

Overall conclusion

The mathematical methods described are incoherent, the supplementary data is not that used to create the figures and equations in the paper, and the supplementary code does not implement the methods described. The paper is methodologically unsound and not reproducible.


Aitchison J. 1986. The statistical analysis of compositional data. London: Chapman and Hall

Aitchison J. 2005. A concise guide to compositional data analysis. University of Glasgow,

Aptroot M. 2016. Yiddish Language and Ashkenazic Jews: A Perspective from Culture, Language, and Literature. Genome Biology and Evolution 8:1948-1949.

Das R, Wexler P, Pirooznia M, Elhaik E. 2016. Localizing Ashkenazic Jews to Primeval Villages in the Ancient Iranian Lands of Ashkenaz. Genome Biology and Evolution 8:1132-1149.

Das R, Wexler P, Pirooznia M, Elhaik E. 2016. Responding to an enquiry concerning the geographic population structure (GPS) approach and the origin of Ashkenazic Jews - a reply to Flegontov et al. arXiv 1608.02038.

Elhaik E, Tatarinova T, Chebotarev D, Piras IS, Maria Calò C, De Montis A, Atzori M, Marini M, Tofanelli S, Francalacci P, Pagani L, Tyler-Smith C, Xue Y, Cucca F, Schurr TG, Gaieski JB, Melendez C, Vilar MG, Owings AC, Gómez R, Fujita R, Santos FR, Comas D, Balanovsky O, Balanovska E, Zalloua P, Soodyall H, Pitchappan R, GaneshPrasad A, Hammer M, Matisoo-Smith L, Wells RS. 2014. Geographic population structure analysis of worldwide human populations infers their biogeographical origins. Nature Communications 5:3513.

Flegontov P, Kassian A, Thomas MG, Fedchenko V, Changmai P, Starostin G. 2016. Pitfalls of the Geographic Population Structure (GPS) Approach Applied to Human Genetic History: A Case Study of Ashkenazi Jews. Genome Biology and Evolution 8:2259-2265.

Pearson K. 1897. Mathematical Contributions to the Theory of Evolution - On a Form of Spurious Correlation Which May Arise When Indices Are Used in the Measurement of Organs. Proceedings of the Royal Society of London 60:489-498.

Friday, 31 January 2014


I've set up this blog as a place to help me gather my thoughts on things that I read or see that relate to my interests. The topics covered are likely to include dating in archaeology and palaeoanthropology, applications of quantitative methods to study the past, historical demography, chemical and isotopic analysis in archaeology with forays into other aspects of archaeology, history and genealogy.

For several years I've been putting off having a blog because I didn't think I could sustain the flow of new writing, but I've finally succumbed. I don't promise to post often or on any regular schedule.