How to know if an individual specimen may be a member of a given reference population?

If the reference population sample is large enough, it might also be possible to apply any (robust) method of outlier detection to the whole data, i.e. the \(n+1\) values composed of the reference sample and the unknown individual. If the unknown individual is detected as an outlier within this sample, this may constitute an evidence against the hypothesis that the unknown individual is a member of the reference sample.

A first possible approach would be to perform a PCA on the whole dataset². For this example, we will use an example from the Goldman Data Set, freely available online [1]. Our reference population sample will be composed of 29 male individuals from the Indian Knoll population sample (US-KT, 5500–3700 BP), and our target individual will be a male individual from Tigara population sample (US-AK, 750–300 BP). Thus, we hope to demonstrate that the target individual is not a member of the reference sample.

Six variables of interest will be considered here (various femoral, humeral and radial bone measurements). Their corresponding short codes in the Goldman Data Set are LFML, LHML, RTML, LHMLD, LFAPD, LFMLD respectively.

library(dplyr)
## Load dataset:
gds <- read.csv("http://web.utk.edu/~auerbach/Goldman.csv",
                dec = ".", na.strings = "",
                stringsAsFactors = TRUE,
                fileEncoding = "macintosh")
## Select the variables of interest and keep only male individuals:
gds <- subset(gds, Sex == "0")
gds <- select(gds, NOTE, LFML, LHML, RTML, LHMLD, LFAPD, LFMLD)
## Select a reference sample and a target individual:
ref <- subset(gds, NOTE == "Indian Knoll")
target <- gds[1, ]
## Merge them in one dataframe:
dat <- rbind(target, ref)
## Drop unused factor levels:
dat <- droplevels(dat)
## Summarize data:
summary(dat)

                      NOTE         LFML            LHML            RTML      
 Indian Knoll           :29   Min.   :404.0   Min.   :292.0   Min.   :337.5  
 Tigara - Point Hope, AK: 1   1st Qu.:424.8   1st Qu.:305.5   1st Qu.:355.5  
                              Median :441.0   Median :312.5   Median :367.0  
                              Mean   :439.9   Mean   :318.0   Mean   :365.8  
                              3rd Qu.:457.8   3rd Qu.:327.5   3rd Qu.:373.0  
                              Max.   :474.0   Max.   :354.0   Max.   :401.0  
                                              NA's   :1       NA's   :3      
     LHMLD           LFAPD           LFMLD      
 Min.   :16.09   Min.   :23.92   Min.   :20.66  
 1st Qu.:18.47   1st Qu.:25.66   1st Qu.:22.55  
 Median :19.10   Median :26.68   Median :23.58  
 Mean   :19.16   Mean   :26.81   Mean   :23.91  
 3rd Qu.:20.29   3rd Qu.:28.09   3rd Qu.:24.59  
 Max.   :24.30   Max.   :31.39   Max.   :29.77

There are only very few missing values here, so that we can simply ignore them and use the default method, i.e. the replacement by the column means³, and perform a PCA:

library(FactoMineR)
## Compute (silently) the PCA:
res_pca <- PCA(dat, quali.sup = 1, graph = FALSE)
## Display a PCA plot:
plot(res_pca, habillage = 1, graph.type = "classic")

As it can be seen on figure 1, the null hypothesis that the target individual (in red) belongs to the reference population sample (in black) is rather unlikely (well, in this toy example, we know that it’s really false!). Adding a 95% data ellipse around the reference population sample is a reasonable way to better visualize this incompatibility:

library(car)
## Display a PCA plot:
plot(res_pca, habillage = 1, graph.type = "classic")
## Extract the principal coordinates of the reference sample
## (we remove the first line, since it corresponds to the target individual):
coor <- res_pca$ind$coor[-1, ]
## Add a 95% data ellipse:
dataEllipse(x = coor[, 1], # PC1
            y = coor[, 2], # PC2
            levels = 0.95, # 95% ellipse
            col = "black", add = TRUE)

The figure 2 is eloquent, and brings an additional proof against the hypothesis that the target individual might be a member of Indian Knoll population sample: it is clearly not within the region of the 95% “most credible” individuals for this population sample. Some remarks, however:

1. This discussion must take into account the percentage of variance captured by the first two principal components (here, it is large enough: it is not mandatory to inspect the following principal components).
2. You can try with other confidence levels in order to get even an stronger evidence!
3. It is not sufficient to say that an individual may not be part of a reference sample: it is important to precise why (i.e., which variables or combination of variables observed on this individual are really unusual). Inspecting the correlation circle of the PCA should be instructive; applying univariate methods (as exposed in the previous section) on well-chosen variables should also help.

As for the univariate case, it might be possible to apply any multivariate outlier detection method to the whole dataset, including both the reference sample and the target individual.

The concept of typicality probability allows to get quantitative information about the hypothesis that an individual belongs to a reference population [3,5], in a way that differs from other methods. This method is particularly useful if there are several population samples as candidates: does the unknown individual belong to population A or population B ?

If the question can really be formulated like this, a common way to solve the problem is to apply a linear discriminant analysis. This will provide posterior probabilities \(p_A, p_B\) for the alternative “the individual comes from population A / population B”. Here, those two hypotheses are mutually exclusive and one of them is necessarily true, so that \(p_A + p_B = 1\).

This framework may be inconvenient in many situations: it makes sense only when there are no other possible reference population for the unknown individual: it is strictly impossible that the individual belong to a third population C. In concrete applications, this cannot always be completely excluded. The goal would be to get a quantitative assessment of the hypotheses “the individual may belong to pop A” / “the individual may belong to pop B”, taking into account that:

• the individual may be incompatible with both populations A and B;
• the measurements of the individual may be perfectly compatible with both A and B (so that none of the two hypotheses can really be discarded).

In other words, the respective credibility of both hypotheses (\(p_A, p_B\)) will not necessarily sum up to 1.

The concept of typicality probability offers an answer to this problem. It might be seen as a multivariate version of the probabilistic distance.

How to compute the typicality probability for a given individual \(x\) and a reference sample? Full details can be found in the seminal article from Wilson [5]. Here is a short summary for the case where an individual is to be compared to one single reference sample:

1. It is assumed that the data within each group follow a multivariate normal distribution.
2. Compute the (squared) Mahalanobis distance between \(x\) and the centroid of the reference sample \(m\): \[ D^2 = (x - m)' V^{-1} (x - m)\] where \(V\) is the covariance matrix of the reference sample.
3. It can be proved that, if \(p\) is the number of variables and \(n\) the number of individuals in the reference sample: \[f_\text{obs} = \frac{n(n-p)}{p(n-1)(n+1)} D^2 \sim \mathcal{F}(p, n-p) \]
4. Thus, the typicality probability is defined as \(1 - F(f_\text{obs})\), where \(F\) is the distribution function of a Fisher distribution with \(p\) and \(n-p\) degrees of freedom.

From its definition, we can see that the typicality probability is the probability to find within the reference population an individual as extreme as, or even more extreme than, the target individual.

Let’s retrieve the example from section 2.1, but we’ll add a second reference population sample: the foragers from Ipituaq (US-AK, 1500–1100 BP). One could expect the target individual to be closer to Ipituaq population than Indian Knoll population, but the target individual might also be far from the two populations.

First, we build a new dataset:

## Add a second population sample from the Goldman Data Set:
dat <- rbind(dat, subset(gds, NOTE == "Ipituaq - Point Hope, AK"))
## Drop unused factor levels, and remove all missing values:
dat <- droplevels(dat)
dat <- na.omit(dat)
## Summarize the data:
summary(dat)

                       NOTE         LFML            LHML            RTML      
 Indian Knoll            :25   Min.   :385.0   Min.   :275.0   Min.   :305.0  
 Tigara - Point Hope, AK : 1   1st Qu.:407.5   1st Qu.:301.0   1st Qu.:337.5  
 Ipituaq - Point Hope, AK:15   Median :429.0   Median :308.0   Median :360.0  
                               Mean   :429.3   Mean   :310.4   Mean   :354.2  
                               3rd Qu.:451.0   3rd Qu.:321.5   3rd Qu.:371.5  
                               Max.   :513.5   Max.   :352.5   Max.   :401.0  
     LHMLD           LFAPD           LFMLD      
 Min.   :16.09   Min.   :23.92   Min.   :20.66  
 1st Qu.:18.52   1st Qu.:26.06   1st Qu.:23.03  
 Median :19.53   Median :27.21   Median :24.64  
 Mean   :20.10   Mean   :27.45   Mean   :25.16  
 3rd Qu.:21.36   3rd Qu.:28.85   3rd Qu.:27.10  
 Max.   :25.91   Max.   :31.71   Max.   :29.96

The R function typprobClass() from the R package Morpho can do what we want: calculate the typicality probability of a given individual x for several groups.

library(Morpho)
## Isolate the target individual, which is in the first row:
indiv <- as.vector(dat[1, -1]) # do not consider the first column here
## Isolate the numeric data that constitute the reference sample,
## i.e. all the dataframe except the first row (target indiv)
## and the first column (grouping factor):
ref <- dat[-1, -1]
## For a greater ease of use below, keep apart the grouping factor:
pop <- dat[-1, "NOTE"]
## Compute typicality probabilities:
typprobClass(x = indiv, data = ref, groups = pop,
             method = "wilson")$probs

  Indian Knoll Ipituaq - Point Hope, AK
1  0.001216739                0.1980937

This suggests that, based on those six bone measurements, the target individual is certainly not a member of Indian Knoll (with a typicality probability as low as 0.001). However, his bone measurements might be compatible with those usually observed within the population of Ipituaq, since the typicality probability associated to this population sample is approximately 0.20.