Linköping University | Department of Computer and Information Science
Master’s thesis, 30 ECTS | Statistics and Machine Learning
2021 | LIU-IDA/STAT-A--21 /018--SE
Bayes Factors for the Proposition of
a Common Source of Amphetamine
Seizures.
Yash Pawar
Supervisor : Anders Nordgaard Examiner : Krzysztof Bartoszek
Upphovsrätt
Detta dokument hålls tillgängligt på Internet - eller dess framtida ersättare - under 25 år från publicer-ingsdatum under förutsättning att inga extraordinära omständigheter uppstår.
Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner, skriva ut enstaka ko-pior för enskilt bruk och att använda det oförändrat för ickekommersiell forskning och för undervis-ning. Överföring av upphovsrätten vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av dokumentet kräver upphovsmannens medgivande. För att garantera äktheten, säker-heten och tillgängligsäker-heten finns lösningar av teknisk och administrativ art.
Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i den omfattning som god sed kräver vid användning av dokumentet på ovan beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan form eller i sådant sammanhang som är kränkande för upphovsman-nens litterära eller konstnärliga anseende eller egenart.
För ytterligare information om Linköping University Electronic Press se förlagets hemsida http://www.ep.liu.se/.
Copyright
The publishers will keep this document online on the Internet - or its possible replacement - for a period of 25 years starting from the date of publication barring exceptional circumstances.
The online availability of the document implies permanent permission for anyone to read, to down-load, or to print out single copies for his/hers own use and to use it unchanged for non-commercial research and educational purpose. Subsequent transfers of copyright cannot revoke this permission. All other uses of the document are conditional upon the consent of the copyright owner. The publisher has taken technical and administrative measures to assure authenticity, security and accessibility.
According to intellectual property law the author has the right to be mentioned when his/her work is accessed as described above and to be protected against infringement.
For additional information about the Linköping University Electronic Press and its procedures for publication and for assurance of document integrity, please refer to its www home page: http://www.ep.liu.se/.
Abstract
This thesis sets out to address the challenges with the comparison of Amphetamine material in determining whether they originate from the same source or different sources using pairwise ratios of peak areas within each chromatogram of material and then mod-eling the difference between the ratios for each comparison as a basis for evaluation. The evaluation of an existing method that uses these ratios to determine the sum of significant differences between each comparison of material that is provided is done. The outcome of this evaluation suggests that there the distributions for comparison of samples originating from the same source and the comparison of samples originating from different sources have an overlap leading to uncertainties in conclusions. In this work, the differences between the ratios of peak areas have been modeled using a feature-based approach. Because the feature space is quite large, Discriminant Analysis methods such as Linear Discriminant Analysis (LDA) and Partial least squares Discriminant Analysis (PLS-DA) have been implemented to perform classification by dimensionality reduction. Another popular method that works on the principle of nearest centroid classifier called as Nearest shrunken centroid is also applied that performs classification on shrunken centroids of the features. The results and analysis of all the methods have been performed to obtain the classification results for classes+1 (samples originate from the same source) and ´1 (samples originate from different sources). Likelihood ratios of each class for each of these methods have also been evaluated using the Empirical Cross-Entropy (ECE) method to determine the robustness of the classifiers. All three models seem to have performed fairly well in terms of classification with LDA being the most robust and reliable with its predictions.
Acknowledgments
I take this opportunity to express my sincere gratitude to all the people who have helped me with their valuable suggestions, expertise, and comments that have been instrumental in writing my thesis.
First, I would like to thank my commissioner for the project and my supervisor at Linköping University Dr. Anders Nordgaard for giving me this opportunity to work on this project. His constant guidance was invaluable in addressing the research questions and methodologies. I would also like to thank my examiner Dr. Krzysztof Bartoszek whose feedback and suggestions have helped me improve my thesis work.
Special thanks to Simon Dunne from the National Forensic Centre for providing me with the necessary resources and the background information required for the project.
Finally, I would like to thank my parents Dilip Pawar and Jyotsna Pawar, who have al-ways been those strong pillars of support that have helped me get through stressful times and have provided me with wise counsel whenever required.
Contents
Abstract iii
Acknowledgments iv
Contents v
List of Figures vi
List of Tables vii
1 Introduction 1
1.1 Background . . . 1
1.2 Aim . . . 2
1.3 Research tasks . . . 2
1.4 Bayes Factor . . . 3
1.5 Validation of a suggested method (Score based evaluation) . . . 5
2 Data 7 2.1 Data Description . . . 7
2.2 Data preprocessing . . . 9
3 Theory and Methods 12 3.1 Linear Discriminant Analysis . . . 12
3.2 Partial Least Squares Discriminant Analysis (PLS-DA) . . . 13
3.3 Nearest Shrunken Centroid . . . 15
4 Results 16 4.1 LDA results: . . . 16
4.2 PLS-DA results . . . 17
4.3 Nearest Shrunken Centroid (NSC) results . . . 20
4.4 Empirical Cross Entropy (ECE) . . . 21
List of Figures
1.1 Distribution of sum of significant differences for linked and unlinked samples . . . 6
2.1 Chromatogram © National Forensic Centre, Linköping, Sweden . . . 7
2.2 Marginal distribution of features . . . 10
2.3 Comparison of Kurtosis . . . 11
4.1 Variance of features . . . 17
4.2 PLS-DA predictions (ncomp=1,5) . . . 18
4.3 PLS-DA predictions (ncomp=10,20) . . . 18
4.4 PLS-DA misclassification rate . . . 19
4.5 PLS-DA explained variance for X (predictors) . . . 19
4.6 PLS-DA explained variance for Y (response) . . . 20
4.7 ECE (Score based Evaluation) . . . 23
4.8 ECE (LDA) . . . 23
4.9 ECE (PLS-DA) . . . 24
List of Tables
1.1 Scale level based on likelihood ratio . . . 4
2.1 Impurities observed in Amphetamine material . . . 8
4.1 Accuracy based on number of features . . . 17
4.2 Amount of shrinkage∆ . . . 21
1
Introduction
1.1
Background
The production and distribution of illicit drugs are a common yet one of the most significant problems prevailing in Europe. The chemical analysis of such materials is very important to unravel a network of drug distribution. At NFC (National Forensic Centre) in Sweden, one such drug, amphetamine, is in its production subjected to a precipitation process where the impurities within the material emerge. These impurities are observed as peak areas on a chromatogram which is observed as an output from Gas Chromatography analysis[16]. A chromatogram is a graphical representation of the concentration of different compounds that are obtained separately as a result of the chromatography of a material. These observed peak areas constitute a kind of chemical fingerprint of the material. The peak areas of impurities can be used to find the similarities between different samples of amphetamine seizures.
Several statistical and non-statistical methods have been used to determine the degree of similarity between different seizures of drugs. Bolck et al.[3] determined the distribution of distances (Pearson correlation distance or Euclidean distance) for linked and unlinked samples of MDMA (3,4-methylenedioxymethamphetamine) tablets. Further, likelihood ratios were calculated based on these distributions. In another approach, Fischer Linear Discriminant Analysis was used by Janhunen et al.[11] to compare street samples of heroin, Neural networks[4] was used for pattern recognition of cocaine samples. The peak areas of impurities are often referred to as the features in comparison to Amphetamine material, However, pre-treatment of peak areas is done using various methods to perform numerically stable calculations. Andersson et al.[1] applied pre-treatment by performing normalization to the sum of peak areas within each sample and further calculated either, the fourth root or logarithm of the normalized data. Power et al.[14] performed weighting and normalization of the data of Amphetamine materials.
At NFC, a new method of using these peak areas has been introduced. It is built on using pairwise ratios of the peak areas of a sample and compares these ratios between two samples suspected to have a common origin. However, a full study of the statistical properties of such ratios and how they should be compared has not yet been made.
1.2. Aim
1.2
Aim
The method that is currently used at NFC suggests comparing the impurity profile (set of peak areas) between two materials for which the question of a common origin has been raised. The conditions in which the chemical analysis is performed are not always the same which leads to differences between the peaks i.e. two materials can have different chemical profiles not because they are obtained from different sources but because of the difference in the conditions of chemical analysis which can lead to incorrect conclusions.
The proposed method suggests the use of pairwise ratios between the peak areas of all the possible combinations of peaks of the impurities within the same chromatogram. It has been observed that the covariance structure of the ratios of peak areas is less sensitive to the differences in conditions of chemical analysis and the ratios constitute much more stable features than the peak areas themselves.
For each material that is analyzed, these ratios are calculated. The comparison between different materials is performed by modeling ’the difference between the ratios’ so that Bayes factors (or likelihood ratios) can be obtained for making inference about the same or different origin of different materials can be made.
1.3
Research tasks
1. Examine the potentiality of using ’the difference between the ratios’ as features to de-termine whether the materials originate from the same or different sources.
2. Recommend a method (or methods) such that it is possible to obtain the Bayes factors (or likelihood ratios) that can be compared with the precision that is defined by the scale of conclusions used at NFC.
3. Comparison of different approaches and investigate the reason for potential different (or similar) behavior.
1.4. Bayes Factor
1.4
Bayes Factor
At NFC, the evaluation of the available forensic evidence is done with the use of likeli-hood ratios. The likelilikeli-hood ratios determine the relative strength of the evidence for the proposition implicitly forwarded by the ’prosecutor’ against the alternative proposition. For the comparison of Amphetamine materials, the proposition (hypothesis) forwarded by the prosecutor shall be referred to as:
Hm: Two amphetamine materials originate from the same batch of manufacturing.
and the alternative proposition (hypothesis) as forwarded by the defense shall be referred to as:
Ha: Two materials originate from different batches of manufacturing.
The hypotheses considered in this case are assumed to be mutually exclusive and exhaus-tive. For a pair of hypothesis Hmand Ha, the prior odds are defined as the ratio of Pr(Hm)
and Pr(Ha)and the posterior odds are defined as the ratio of Pr(Hm|D)and Pr(Ha|D)where
D is the available data which can specifically be referred to as the impurity profile of the observed amphetamine material for which the comparison is made. The Bayes factor B can then be defined the ratio of posterior odds to the prior odds.
B= Pr(Hm|D)/Pr(Ha|D)
Pr(Hm)/Pr(Ha) (1.1)
When both the hypotheses are simple, the Bayes factor B simplifies to:
V= L(D|Hm)
L(D|Ha) (1.2)
Where, L(D|HM)and L(D|HM)are the likelihoods for each of the hypothesis Hmand Ha
respectively, considered to be true.
Scale of conclusions
In order to support the propositions forwarded by either the ’Prosecutor’ or the ’Defense’, the strength of evidence has to be determined. The evidence is supported by probabilis-tic reasoning. However, it can be difficult to interpret the probabilities for the end users outside the field of expertise. At NFC, an ordinal scale of conclusions is used to report the strength of evidence. The scale of conclusions is based on the value of likelihood ratios (or Bayes Factor) obtained. The levels correspond with intervals of the likelihood ratio V [13]. The scale of conclusions with their corresponding significance used by The Swedish National Laboratory of Forensic Science (SKL)[17] are given below:
1. Level +4: The results are at least 1,000,000 times more probable if the main hypothesis is true compared to if the alternative hypothesis is true.
2. Level +3: The results are at least 6000 times more probable if the main hypothesis is true compared to if the alternative hypothesis is true.
1.4. Bayes Factor
Table 1.1: Scale level based on likelihood ratio Interval Scale Level
V ď 10´6 -4 10´6ďV ď 1/6000 -3 1/6000 ď V ď 1/100 -2 1/100 ď V ď 1/6 -1 1/6 ď V ď 6 0 6 ď V ď 100 +1 100 ď V ď 6000 +2 6000 ď V ď 106 +3 106ďV +4
3. Level +2: The results are at least 100 times more probable if the main hypothesis is true compared to if the alternative hypothesis is true.
4. Level +1: The results are at least 6 times more probable if the main hypothesis is true compared to if the alternative hypothesis is true.
5. Level 0: The results are equally probable if the main hypothesis is true compared to if the alternative hypothesis is true.
6. Level -1: The results are at least 6 times more probable if the alternative hypothesis is true compared to if the main hypothesis is true.
7. Level -2: The results are at least 100 times more probable if the alternative hypothesis is true compared to if the main hypothesis is true.
8. Level -3: The results are at least 6 000 times more probable if the alternative hypothesis is true compared to if the main hypothesis is true.
9. Level -4: The results are at least 1 000 000 times more probable if the alternative hy-pothesis is true compared to if the main hyhy-pothesis is true.
Note: Here, the main hypothesis refers to the hypothesis forwarded by the prosecutor i.e. Hmand the alternative hypothesis refers to the hypothesis refers to the hypothesis forwarded
1.5. Validation of a suggested method (Score based evaluation)
1.5
Validation of a suggested method (Score based evaluation)
A comprehensive study regarding the method of using the pairwise ratio of peak areas that are currently adopted at NFC is done to determine whether the results obtained can be appreciated based on the scale of conclusions that are currently used at NFC to determine the strength of evidence.
Significant Differences
As described in the previous sections, the data that is used for experimentation in this thesis work is based on seizures of Amphetamine material. However, this study was performed on the data that is observed as a result of chemical profiling of seizures of Cannabis hash material. As a result of chemical profiling based on the chromatogram, the peak areas of 17 impurities that are present in each sample were provided.
The pairwise ratios of peak areas that are present in a chromatogram from within each sample were calculated and a matrix of ratios (also referred to as quotients) is constructed (For 17 peak areas, 17 rows, and 17 columns). Comparisons between different samples along with their corresponding level of likelihood ratios are available as anticipated ground truth. For each comparison between different samples that were specified, the absolute difference between the corresponding quotient matrices was calculated resulting in a difference matrix. From the difference matrix, the values corresponding to a significant difference were determined based on the comparison with the error compensated values for each ratio that was obtained. The error compensated values were determined for each ratio by taking into account the uncertainties in the measurements due to various conditions such as the amount of time for which the sample was stored, the temperature at the time of measurement, etc. The Sum of all the significant differences (also regarded as scores) that were present in the difference matrix was calculated and calibrated on a scale according to the level of likelihood ratios that is available (as explained above). The newly calibrated scale based on the scores was then used to determine the extent of similarity or dissimilarity between the samples that were compared[22, 2].
Evaluation
To evaluate the discriminatory power of the obtained calibrated scale based on the scores, the available comparisons between different samples along with their scores and the corre-sponding scale level based on the intervals in which their likelihood ratios lie are taken into consideration. The scores corresponding to the positive values (+1 to +4) and the scores corresponding to the negative values (´1 to ´4) on the scale level are categorized as the samples deemed to originate from the same source and the samples deemed to originate from different sources respectively. The scores corresponding to the level 0 are randomly categorized into two groups. The distribution of the scores corresponding to each of these two groups is determined using a kernel density estimate (KDE).
Figure 1.1 shows the plot of distributions of the sum of significant differences for the samples that are deemed to originate from the same manufacturing batch and the samples that are deemed to originate from different manufacturing batches. The overlapping region between the two distributions resulted in uncertainty in conclusions.
1.5. Validation of a suggested method (Score based evaluation)
2
Data
2.1
Data Description
For each material that is seized and extracted, impurities are observed as peak areas in a chromatogram. The chromatogram is observed as a result of a process called Gas chromatog-raphy combined with Mass selective Detector (GC-MSD) or with Flame Ionization Detector (GC-FID) [16]. In a chromatogram from an analysis of an amphetamine material, peaks of a set of 30 monitored impurities (chemical compounds) are observed as shown in Figure 2.1.
Table 2.1 shows the chemical compounds that are observed as impurities in Am-phetamine Material:
2.1. Data Description
Table 2.1: Impurities observed in Amphetamine material
No. Chemical compounds
1 Ketoxime 1 2 Ketoxime 2 3 4-Methyl-5-phenylpyrimidine 4 Unknown C 5 N-Benzylpyrimidine 6 N-Acetylamphetamine 7 N-Formylamphetamine 8 1,2-Diphenyletylamine 9 N,N-Dibenzylamine 10 1,2-Diphenylethanone 11 Benzylamphetamine 12 DPPA 13 DPIA 1 14 DPIA 2 15 alfa-Methyldiphenetyletylamine 16 DPIMA 1 17 DPIMA 2 18 Unknown A2 19 Naphthalene 1 20 Unknown A3 21 Naphthalene 2 22 N-Benzoylamphetamine 23 Unknown B2 24 2-Oxo 25 2,6-Dimethyl-3,5-diphenylpyridine 26 2,4-Dimetyl-3,5-diphenylpyridine 27 Pyridine 7 and 14 28 2,6-Diphenyl-3,4-dimethylpyridine 29 DPIF 1 30 DPIF 2
The peak areas of each of the chemical compounds (impurities) listed in table 2.1 that are observed in a chromatogram for each material are regarded as features within the data. A set of these peak areas are regarded as the impurity profile of each sample. The data consists of impurity profiles of 741 samples.
A comparison is performed between samples to determine whether they have same ori-gins or different oriori-gins. Various evaluation methods are used for the comparison of samples by considering a pair of samples at a time.
2.2. Data preprocessing
2.2
Data preprocessing
The comparisons between different samples that are available in the data correspond to levels +4,+3,-4. The levels +4 and +3 strongly support the hypothesis stating both the samples are linked (i.e. originate from the same manufacturing batch)and the level -4 strongly supports the hypothesis that states both the samples are unlinked (i.e. originate from different man-ufacturing batches) . As a first step, the given levels are translated to ’+1’ for the levels +4 and +3 and ’-1’ for level -4. After the data has been transformed to new classes (i.e. +1 and -1), the data consists of 499 samples that correspond to class +1 (35.61% of total samples) and 902 samples that correspond to class -1 (64.39% of total samples). When the data is split into training and test sets, the same ratio has been maintained.
The calculation of pairwise ratios of peak areas for all possible combinations of reference peaks within the same chromatogram is performed using the formula:
Kij=
Pi
Pi+Pj (2.1)
where, Piand Pjare the peak areas within the same chromatogram. The k values obtained
are self normalizing in nature and lie between 0 and 1.
As a result, for each sample that represents an impurity profile of a material, a total of 900 ratios are obtained. Ratios for 741 samples for which the peak areas are available, are calculated. For the 1401 samples in comparison, the difference between the corresponding ratios is calculated.
Thus, a sparse matrix comprising of 1401 rows (comparison between samples) and 900 columns (ratios) is obtained. The matrix will be referred to as di f fkmatrix henceforth in the
report . The variance for each feature is calculated which resulted in 44 features having zero variance (values within the feature are constant), Such features are eliminated from the data.
Distribution of features
The marginal distribution of features that are obtained suggest a heavy-tailed distribution centered around zero. Figure 2.2 describes the marginal distribution of few of the feature variables.
A measure describing the nature of distribution whether it is heavy-tailed or light-tailed as compared to a normal distribution is kurtosis. The kurtosis for a standard normal distri-bution is 3. A heavy tailed distridistri-bution has a high kurtosis whereas the value for a light tailed distribution is low. Thus, the excess kurtosis gives a measure in comparison to the kurtosis of a standard normal distribution. The excess kurtosis is defined as:
γ2=m4/m22´3 (2.2) where, m4= n ÿ i=1 (xi´¯x)4/n (2.3) and, m2= n ÿ i=1 (xi´¯x)2/n (2.4)
2.2. Data preprocessing
Figure 2.2: Marginal distribution of features
¯x is the mean and n is the sample size.
Lambert W transformation:
Very many statistical methods have been developed under the assumption that the data is Gaussian. In practice, however, that is not necessarily the case with the data that is observed in nature as it can possess some kind of asymmetry or skewness. Thus, inferences made from the models based on this type of data can lead to inaccuracies.
According to Goerg et.al. in [8], a generalized LambertW ˆ Fxrandom variable Z can be
defined for a continuous random variable X for which the Cumulative Distribution Function (CDF) is defined as Fxand is governed by a tail parameter δ as shown in equation 2.5 where
U = (X´µX
σX )and X is continuous random variable such that X Fx(x|β)and µX and σX are
the mean and the standard deviation of the random variable X.
For the values δ ě 0 the equation 2.5 becomes a bijective function. When X ãÑN (µX, σX2)
the transformation becomes Lambert W ˆ Gaussian and δ=0 results in Z being a Gaussian variable. The increase in value of δ, leads to heavier tails. The inverse transformation of the equation 2.5 is described in section 2.3 of [7] which gives a function to gaussianize a
2.2. Data preprocessing
The methods used for experimentation in this thesis work assumes data to be Normally distributed. Thus, a transformation Lambert W is used to transform the distribution of ob-served marginal features in the data to the distribution for which the γ2value as described
in equation 2.2 is close to 0 (which resembles the kurtosis of standard normal distribution). This, however, does not guarantee that the distribution is Gaussian but can work as a base-line for a Gaussian sample[7]. The implementation has been done in R using the function
1gaussianize1 from the package1Lambert W1[9]. Figure 2.3 shows the comparison of excess
kurtosis of the original variables and the Gaussianized variables.
3
Theory and Methods
3.1
Linear Discriminant Analysis
Linear Discriminant Analysis (LDA) is one of the most successful and popular methods when it comes to classification for a high dimensional data. LDA has been developed for explanatory variables that follow normal distributions. Thus, results of LDA are reliable when the assumption of normality in the data is fulfilled. Another assumption that LDA makes is that the co-variance matrices for all the classes are equal i.e.
Σ1=Σ2=...=Σk (3.1)
LDA identifies the linear features that maximizes the between-class separation while minimizing the within class separation. Consider a training set x1, ..., xn consisting of N samples. Each sample xiis a column vector of length d. The training samples belong to two
classes C1 and C2. The sample points are projected on a line that maximize the separation
between two classes C1and C2.
The projection of points are given as:
ai=βTXi (3.2)
3.2. Partial Least Squares Discriminant Analysis (PLS-DA)
The idea is to find a projection for which the class means are as far apart as possible while minimizing the within class variances [18]. This can be achieved by:
max
β:||β||=1
(µ1´ µ2)2
s21+s22 (3.5)
Implementation
The method has been implemented on the Gaussianized features from the di f fk matrix.
There are a few features in the data that have very low variance and appear to be constant within the data usually tend to have a high co-linearity amongst them. Thus, appropriate number of features are selected based on the marginal variance as a threshold.
The variance is used as threshold by performing a grid search over all the set of values present in the grid. A sequence of values are used to perform a grid search ranging from Var=0.001 to Var=0.4. LDA is performed based on the features that have variance higher than the value present in the grid. The evaluation is performed for the model using a leave one out cross validation method. The features corresponding to the model with the highest accuracy (or the least misclassification error) are used to train the model.
For training and testing of the LDA model, the data is split into 80% and 20% respectively. The LDA model is fitted with the training data where the prior probabilities of both the classes
+1 and ´1 are set to Pr(Hm) =Pr(Ha) =0.5 making it a uniform prior. Based on the LDA
model that is fitted, the classes are predicted from the new data and the posterior probabilities for each sample are obtained. All the calculations and experimentation is done in R.
3.2
Partial Least Squares Discriminant Analysis (PLS-DA)
The method PLS-DA is a variant of Partial Least Squares Regression (PLS-R) when the re-sponse variable is categorical. PLS-R is a projection method that aims towards finding, for an M-dimensional X variable (consisting of N samples) , a new A-dimensional hyperplane in such a way that the correlation between these new projections and the response variable Y is maximized[23].
The method is widely used in situations where the data consists of large number of features in comparison with the number of samples and are possibly, highly correlated. This situation is commonly observed when dealing with Chemometric analysis (a branch in chemistry that employs statistical, mathematical methods to derive conclusions)[10].
The method is used to find the covariance structure between the predictor and the re-sponse variables by projecting them onto a new space. In short, PLS-DA performs a linear transformation that projects the original data into a lower dimensional space. The reduced dimensions in this case are referred to as principal components (PCs). The PCs are formed in such a way that covariance between the original data and the response variable is preserved with the first PC explaining the highest amount of covariance [6].
3.2. Partial Least Squares Discriminant Analysis (PLS-DA)
The decomposition of both the predictor variable X and the response variable Y is given by the equations:
X=TPT+E (3.6)
Y=TQT+F (3.7)
where, T is defined as a score matrix and P and Q are the matrices of the loadings of X and Y, respectively. E and F are the error matrices associated with X and Y matrices, respectively.
If a weight matrix W is defined for the decomposition of data matrix X, the score matrix T is given as:
T=XW(PTW)´1 (3.8)
Analogous to linear regression, Partial least square regression can be written as: ˆ
Y=X ˆβ (3.9)
The co-efficient matrix ˆβ can be obtained by replacing the T from equation 3.7 in the equation 3.6. The matrices W, T, P and Q can be obtained using an algorithm nonlinear iterative partial least squares (NIPALS) as stated by Tang et. al.[19].
Implementation
Similar to the implementation of LDA as described in the previous section, Gaussianized features are used to perform Partial least squares discriminant analysis. The PLS-DA model is trained on the training data consisting of 1120 samples and for validation, a leave-one-out cross validation is performed. A linear projection for 20 components is obtained. The testing is done on test data consisting of 281 samples. The classification accuracy for 20 components is retrieved. All the calculations and simulations are done in R and the package mdatools is used for implementation of the PLS-DA model.
3.3. Nearest Shrunken Centroid
3.3
Nearest Shrunken Centroid
The method is primarily devised for profiling of gene expression in a way so as to perform class predictions based on the nearest centroid classifier. Considering the high amount of features in the data, this method shrinks the within class centroid for each gene expression (feature) towards the overall centroid of the data. A threshold is used to determine the shrinkage within each class for each feature. If the distance of the centroid is less than the threshold, it shrinks the feature to zero, thus reducing the effect of noisy features. The threshold is decided by cross-validation[20].
Amount of shrinkage
The amount of shrinkage is determined by performing n-fold cross validation with different values of shrinkage parameter∆ as shown in equation 3.10.
d1ik=sign(dik)(|dik´∆|)+ (3.10)
where, dikis the standardized difference between the centroid for gene expression i within
class k and the overall centroid of gene expression i. The standardization performed in this case is analogous to the one performed in LDA. The classification is done based on the comparison with each of these class centroid. The class whose centroid is the sample closest to, is the predicted class for the sample[20].
Implementation
The implementation of nearest shrunken centroid method is done in R and the package used for implementation is1pamr1[20]. A 10 fold cross validation is performed to determine the
threshold for shrinking the class centroid towards the overall centroid of the data. The train-ing data consisttrain-ing of 1120 samples is used for selection of threshold and fitttrain-ing the model. The prior probabilities of both the classes+1 and ´1 are set to Pr(Hm) =Pr(Ha) =0.5
mak-ing it a uniform prior. The prediction is done for the test data consistmak-ing of 281 samples which resulted in the predicted classes as well as the posterior probabilities for each sample.
4
Results
The results presented in this section are based on the experimentation that is presented in the theory and method section. The experimentation involves implementation of Linear Discriminant Analysis (LDA) and Partial Least Squares Discriminant Analysis (PLS-DA) to determine whether the comparison of samples in question originate from the same manufac-turing batch or from different manufacmanufac-turing batches.
As discussed in the previous sections, the hypothesis for samples originating from same manufacturing batch Hmand the hypothesis for samples originating from different
manufac-turing batches Haare denoted as levels+1 and ´1 respectively for the comparisons between
the samples that are available in the data as the anticipated ground truth. Thus, results of the experimentation are primarily discussed in terms of misclassification error.
In order to obtain the likelihood ratios (as described in equation 1.2), the posterior distri-bution for each of the hypothesis Hmand Hais also calculated. By entering prior odds equal
to 1, instead of taking the prior odds from the ratio of the number of same source cases to the number of different source cases in data. Then the denominator of equation 1.1 will be 1 and the posterior odds equal to the Bayes factor.
The results for the experimentation performed are discussed in the sections that follow.
4.2. PLS-DA results
Figure 4.1: Variance of features
Table 4.1: Accuracy based on number of features Features Accuracy Variance
608 0.8979300 0.001 369 0.9329051 0.045 260 0.9072091 0.089 170 0.8857959 0.133 91 0.8379729 0.177 55 0.7765882 0.221 24 0.7480371 0.265 6 0.7173448 0.309 6 0.7173448 0.353 3 0.5695931 0.397
4.2
PLS-DA results
In this experimentation, the results are obtained by performing a leave-one-out cross valida-tion on the training data and the predicvalida-tion of classes is also performed on the test data. The training and test data are divided into 80-20%. The parameters to be tuned in this method are the number of components. The projections for as many as 20 components are obtained to validate the accuracy of the class predictions and also to examine the variance explained by the components.
Class Predictions based on number of components
The performance of PLS-DA model is evident from the prediction results as shown in figure 4.2 and 4.3. It can be seen that when ncomp=1, the predictions are not accurate. However, with the increase in number of components to ncomp= 5, the accuracy is improved signif-icantly. As the number of components keep on increasing the performance keeps getting better. It is worth noting that for ncomp=10 to ncomp=20 there is no significant improve-ment in the results. Analogously, the decrease in the rate of misclassification is evident with
4.2. PLS-DA results
the increase in ncomp from the figure 4.4.
4.2. PLS-DA results
Figure 4.4: PLS-DA misclassification rate
for the first few components is also poor. However, when we look at the cumulative variance for ncomp = 10 and onward, the explained variance is quite high. In PLS, the components are formed in such a way that the covariance between the linear projections tj = XWj and
response variable Y is maximized. Therefore, the optimal number of j are to be defined so as to ensure reliable predictions. It is important to determine the predictive significance of each added component and to ensure stopping when the number of components start becoming insignificant as it might lead to over-fitting. The components start becoming insignificant after ncomp= 10 as the increase in the amount of variance explained by the following com-ponents is insignificant.
4.3. Nearest Shrunken Centroid (NSC) results
Figure 4.6: PLS-DA explained variance for Y (response)
4.3
Nearest Shrunken Centroid (NSC) results
The class centroid for each gene expression (feature) is shrunk towards the overall centroid of the data by a factor∆ (equation 3.10) using soft thresholding as described in section 3.3. The results for the thresholding that is determined by cross validation are shown in the table 4.2. The classification error is minimum for∆=3.7991 and the number of gene expressions (fea-tures) that appear to be stable for this threshold are 76. Thus, only a 10% of gene expressions (features) contribute towards classification.
4.4. Empirical Cross Entropy (ECE)
Table 4.2: Amount of shrinkage∆ Size Threshold Error
750 0.0000000 0.2383929 590 0.2532747 0.2330357 529 0.5065493 0.2339286 465 0.7598240 0.2357143 404 1.0130987 0.2375000 360 1.2663734 0.2410714 337 1.5196480 0.2419643 292 1.7729227 0.2455357 263 2.0261974 0.2464286 243 2.2794721 0.2446429 217 2.5327467 0.2401786 185 2.7860214 0.2383929 159 3.0392961 0.2348214 132 3.2925708 0.2241071 108 3.5458454 0.2250000 76 3.7991201 0.2196429 65 4.0523948 0.2205357 53 4.3056694 0.2883929 42 4.5589441 0.4080357 27 4.8122188 0.4455357 24 5.0654935 0.4500000 17 5.3187681 0.4357143 12 5.5720428 0.4035714 9 5.8253175 0.3803571 4 6.0785922 0.3642857 4 6.3318668 0.3580357 4 6.5851415 0.3553571 4 6.8384162 0.3553571 2 7.0916909 0.3553571 0 7.3449655 0.3553571
4.4
Empirical Cross Entropy (ECE)
The Empirical Cross Entropy (ECE) is a concept from the field of Information theory that is used to measure the accuracy of the calculated likelihood ratios. The accuracy is measured in terms of average information loss. According to Information theory, the average information obtained as a result of an inferential process is determined by the reduction of loss in the entropy of the system. The entropy gives the measure of uncertainty (or randomness) within a given system. The uncertainty of an event can be determined using the probability of the occurrence of the event. A higher probability signifies the certainty of the occurrence of the event and thus less randomness specified by a low entropy[15].
Based on the obtained likelihood ratios, the posterior probability for the hypothesis Hm
provided the evidence is given as shown in equation 4.1.
P(Hm|D) = LR.P(Hm) P(Ha) 1+LR.P(Hm) P(Ha) (4.1)
4.4. Empirical Cross Entropy (ECE)
where, P(Hm)and P(Ha)are the prior probabilities for hypothesis Hmand Haand LR is
the the likelihood ratio of the hypothesis Hm to the hypothesis Ha given the evidence D as
given in equation 1.2.
Assuming that the prior probabilities P(Hm)and P(Ha)are known, the empirical cross
entropy (equation 4.2) in its general form measures the amount of information needed to know which of the two hypothesis Hmand Hais true.
ECE=´P(Hm) Nm . ÿ jeDm log2P(Hm|Dj)´ P(Ha) Na . ÿ jeDa log2P(Ha|Dj) (4.2) Where,
• Nm: Number of comparison of samples originating from the same origin.
• Na: Number of comparison of samples originating from different origins.
• Dm: Comparison of samples originating from the same origin.
• Da: Comparison of samples originating from different origins.
The ECE plots for the score based evaluation as well as for the LDA method are shown in figure 4.7 and 4.8. The ECE values are computed and plotted as a function of the range of log odds of prior probability θ=P(Hm). The values of prior probability θ are considered for the
range[0, 1]. For each of these prior probabilities, the posterior probability P(Hm|D), and the
likelihood ratio values, the ECE values are computed using equation 4.2 and are represented as a function of P(Hm).
4.4. Empirical Cross Entropy (ECE)
Figure 4.7: ECE (Score based Evaluation)
Figure 4.8: ECE (LDA)
Figure 4.7 shows the ECE plot for the score based evaluation. The red curve represents the ECE values computed for the likelihood ratios of the comparison of samples originating from the same source and the comparison of samples originating from different sources. The higher is the peak of the curve, the more is the uncertainty within the system and therefore, more information is needed to know the true hypothesis[15].
The black curve represents the ECE of the system when the LR (likelihood ratio) is iden-tical to 1 suggesting a neutral system. The ECE value of such a system is simply the cross entropy of the prior probabilities of the hypotheses. If, for a system, the red curve is higher than the black curve, there is a more loss of information when the decisions are based on the
4.4. Empirical Cross Entropy (ECE)
4.5. Performance Comparison
likelihood ratios calculated for such a system than when the decisions have been made based on the prior information. The blue curve represents the ECE for a calibrated system that is based on the Pooled Adjacent Violators Algorithm (PAVA) (see [12] for details).
From the comparison of plots 4.7 and 4.8, it can be seen that the ECE plot for the likeli-hood ratios obtained using LDA method is almost flat suggesting almost a perfect system which requires very little to no additional information to make inference. The discrimination achieved in this system can be assumed to be more reliable and robust as compared to the score based evaluation.
The ECE plots for the likelihood ratios obtained for the comparisons of samples originat-ing from same source and the samples originatoriginat-ing from different sources from PLS-DA and Nearest shrunken centroid model are shown in figure 4.9 and 4.10.
4.5
Performance Comparison
Comparison of results obtained for models LDA, PLS-DA and the Nearest centroid classifier (NSC) are shown in the table 4.3. The performance metrics used for evaluation are accuracy, sensitivity, specificity, precision and the balanced accuracy. All the performance metrics are expressed in percent.
Table 4.3: Classification Results
Performance Metrics LDA PLS-DA NSC
Accuracy 92.53 86.83 74.02
Sensitivity 92.71 94.79 94.79 Specificity 92.43 82.70 63.24
Precision 86.41 73.98 57.23
5
Discussion
This chapter discusses the results and relevance of the methods used to approach the thesis objectives and also the challenges and limitations identified in the implementation of these methods.
5.1
Results
Performance Evaluation
The results of both the discriminant analysis methods (LDA and PLS-DA) for classification have been satisfactory. The LDA model has been optimized by considering the features with variance only above a certain threshold which is obtained through cross-validation over a set of variance values. However, only the marginal variance of features has been taken into consideration. This step has been implemented to eliminate features that appear to be almost constant. Such features lead to a very high colinearity amongst themselves. PLS-DA, on the other hand, tries to find the projections in a low dimensional subspace such that its covariance with the response variable is maximized. As discussed in section 4.2, in PLS-DA implementation, the classification results for a maximum of ncomp = 20 are obtained, of which, ncomp= 10 are selected for the performance evaluation. Both the methods perform classification by performing dimensionality reduction as a first step.
For both the methods, the performance evaluation as shown in table 4.3 is done on the test data set. The accuracy of LDA, however, is better than the accuracy of PLS-DA. Specificity
5.1. Results
NSC is also implemented to perform classification after performing feature reduction by shrinking of centroids. As discussed in section 4.3, only 10% of total features appear to be stable after shrinking the centroids and are used to perform classification. The accuracy obtained for this method is lower than both LDA and PLS-DA. However, the sensitivity is higher than both LDA and PLS-DA. The specificity is very low at 63.24% which means that it would lead to a very high FPR.
Contribution of features
The feature space that is considered in all three methods is quite large. Hence, it is important to observe the features that are common amongst all the methods. The information about the contribution of features amongst different models is quite significant to the experts at NFC as it provides information about the peak areas of the chemical impurities that contribute more towards classification than others. In the case of LDA and PLS-DA, 200 most important features are identified based on the coefficients (as described in equations 3.2 and 3.9) that form their projections, out of which, only 44 of them appear to be the same. In the case of NSC (as shown in table 4.3) 76 gene expressions (features) are identified to be the most stable and have been used to train the model, considering these features to be the most contributing ones. When common features are determined amongst all the three models, only 8 appear to be the same.
Likelihood Ratios
This thesis work is one of the first attempts to use a feature-based evaluation method for modeling the difference between the pairwise ratios peaks areas within a chromatogram for the comparisons between the samples that are available. Because for evaluation of the strength of evidence at NFC, the likelihood ratios are used, it is important to suggest methods that can yield the results in such a way that they can be converted to likelihood ratios.
Each of the methods that are suggested in this thesis work yield the results in the form of a conditional posterior distribution for each class that is specified in the model. The ratio of posterior probabilities for one class to the other (+1 to ´1) for every sample is identical to likelihood ratios given that the prior probabilities for each class are the same.
The ECE plots as shown in section 4.4 describe the accuracy and robustness of likelihood ratios obtained using each of these methods. For LDA, the ECE plot for the obtained likeli-hood ratios appears to be quite stable and can be interpreted to be a near-perfect system as the red curve as shown in figure 4.8 appears to resemble the calibrated system. In the case of PLS-DA, the peak for the red curve as shown in figure 4.9 is higher than LDA but not higher than a neutral system that is shown with the black curve. When the ECE plots likelihood ratios obtained from the NSC method are plotted as shown in figure 4.10, they appear to be higher than the ECE plot for LDA but lower than the ECE plot for PLS-DA. In general, all the three systems can be accounted to be performing fairly well.
The ECE plots can also be used to interpret the discriminatory power of the systems in discussion. Lower is the curve, better is the discriminatory power. In our case, the dis-criminatory power can be defined to be the separation of likelihood ratios for each of the hypotheses. For instance, the likelihood ratio for the samples originating from the same batch (i.e. class+1) is much higher in the case of LDA than the likelihood ratios obtained for samples originating from the same batch in the case of PLS-DA. Similarly, it is much lower for samples originating from different batches (class ´1) in the case of LDA as compared to samples originating from different batches in the case of PLS-DA. Thus, the separation of
5.2. Method
likelihood ratios is much higher for classes+1 and ´1 in LDA as compared to PLS-DA. It can however happen that the accuracy of a model (PLS-DA) is higher than the other (NSC) but the discriminatory power of the other (NSC) is better for the former (PLS-DA).
The results of the ECE plot can also be compared with the score-based evaluation method as described in section 1.4. It can be seen that the separation of likelihood ratios for the score-based method is lower than that of LDA and NSC and is quite comparable to the PLS-DA method. Hence, it can be said that calibrating the scale levels based on the significant scores may lead to inconsistencies in decision-making.
5.2
Method
Gaussianizing the variables as a Pre-processing method for data
There are various reasons for the pre-treatment of the predictor variables present in the data. In our case, the pre-treatment is done to address the assumption of the proposed models that the predictor variables follow Gaussian distributions. The transformation Lambert W enables using the properties of Gaussian distribution for the transformed data for the models that have been developed assuming the distribution of the predictor variables to be Gaussian[7].
Implementation
1. LDA has proved to be one of the most sophisticated classification methods for high-dimensional data. The separation of likelihood ratios achieved for this model is quite large thus leading to a very high discriminatory power. The results however cannot be relied on entirely as it is the performance of a model in terms of classification as well as the separation of likelihood ratios is yet to be observed for the data corresponding to scale levels+2 to ´3 for which the support of the samples originating from same or different batches is less probable as compared to the scale levels provided in the data. (+3,+4, ´4). LDA also faces a problem of collinearity amongst the predictor variables when the number of observations (n) in the data is less than the number of predictor variables (p).
2. PLS-DA often seems to perform better in situations where p ą n. The method also tries to find projections that have higher covariance with the response variable. Although the classification results obtained with projections for ncomps = 10 or more seem to be quite good, the separation between the likelihood ratios as discussed in section 5.1 is very small which might lead to a lack of confidence for the decision-maker to reach conclusions.
5.3. Future Work
or PLS-DA.
5.3
Future Work
There are methods through which the work done in this thesis can further be improved especially considering the feature-based modeling approach for the pairwise ratio of peak areas. The discriminant analysis methods used in this thesis mostly are implemented as a single-stage process. However, considering the high dimensionality of the data, a two-stage process can be used such as Principal Component Analysis-Linear Discriminant Analysis (PCA-LDA) where PCA can be used as an initial step to perform dimensionality reduction and the output scores of PCA can then be used as an input to LDA [21]. Another variant of such implementation is PLS-LDA[19]. Such implementations can prove very powerful for some specific applications especially when p ąą n.
Although the methods suggested in this thesis work provide satisfactory results in terms of classification, the performance of these methods for data with neutral scale levels is yet to have experimented with. Also, the discriminatory power of these models is not yet tested when the raw scale levels ´4 to+4 are used directly instead of the+1 (’Same source’) and ´1 (different sources) with the data. This limitation occurs when the number of observations in the data is insufficient as compared to the number of classes (as in this case it is 9).
6
Conclusion
This thesis work addresses the main challenges of evaluation of the existing score-based method by evaluating the likelihood ratios for the comparison of samples and analyzing them with the help of ECE plots. The outcome of this evaluation shows that the two distribu-tions (Samples originating from same source and the samples originating from different sources) have a significant amount of overlap (figure 1.1) suggesting that there could be uncertainties in classification.
As an alternative to the existing score-based approach, a feature-based modeling ap-proach has been proposed for classification. LDA, PLS-DA, and NSC have been proposed for classification considering the high dimensional feature space. All the methods seem to have performed well with LDA providing the best overall accuracy in terms of classification. PLS-DA and NSC also proved to be good in terms of accuracy with PLS-DA performing slightly better than NSC.
The comparisons between the suggested methods are also done based on the likelihood ratios for each comparison of samples (pairwise). As a result, the extent of separation between likelihood ratios of the two classes (+1 and ´1) is compared to determine the discriminatory power of the methods. LDA has proved to be the most robust of the three in terms of its discriminatory power suggesting very high confidence in the decision. NSC also provides good confidence in the decision based on its discriminatory power although the classification accuracy is not better than LDA and PLS-DA.
Bibliography
[1] Kjell Andersson, Eric Lock, Kaisa Jalava, Henk Huizer, Sten Jonson, Elisabet Kaa, Alvaro Lopes, Anneke Poortman-van der Meer, Erkki Sippola, Laurence Dujourdy, and Johan Dahlén. “Development of a harmonised method for the profiling of am-phetamines VI: Evaluation of methods for comparison of amphetamine”. In: Forensic science international 169.1 (), pp. 86–99.
[2] Anoek Biemans. “Improvements in the Amphetamine Profiling Method at NFC : Semi-automatic chromatogram analysis, reference substances and drying method”. In: (2020).
[3] Annabel Bolck, Céline Weyermann, Laurence Dujourdy, Pierre Esseiva, and Jorrit van den Berg. “Different likelihood ratio approaches to evaluate the strength of evidence of MDMA tablet comparisons”. In: Forensic Science International 191.1 (), pp. 42–51. [4] John Casale and J Watterson. “A Computerized Neural Network Method for
Pat-tern Recognition of Cocaine Signatures”. In: Journal of forensic sciences 38 (Apr. 1993), pp. 292–301.
[5] Byeong Yeob Choi, Eric Bair, and Jae Won Lee. “Nearest shrunken centroids via alter-native genewise shrinkages”. In: PLOS ONE 12.2 (Feb. 2017), pp. 1–18.DOI: 10.1371/
journal.pone.0171068.
[6] Mario Fordellone, Andrea Bellincontro, and Fabio Mencarelli. Partial least squares dis-criminant analysis: A dimensionality reduction method to classify hyperspectral data. 2018. arXiv: 1806.09347 [stat.AP].
[7] Georg Goerg. “The Lambert Way to Gaussianize Heavy-Tailed Data with the Inverse of Tukey’s h Transformation as a Special Case”. In: The Scientific World Journal 2015 (Oct. 2010).
[8] Georg M. Goerg. “Lambert W random variables—a new family of generalizedskewed distributions with applications to risk estimation”. In: The Annals of Applied Statistics 5.3 (2011), pp. 2197–2230.
[9] Georg M. Goerg. LambertW: Probabilistic Models to Analyze and Gaussianize Heavy-Tailed, Skewed Data. R package version 0.6.6. 2021.
Bibliography
[10] KÁROLY HÉBERGER. “Chapter 7 - Chemoinformatics—multivariate mathemati-cal–statistical methods for data evaluation”. In: Medical Applications of Mass Spectrome-try. Ed. by Károly Vékey, András Telekes, and Akos Vertes. Amsterdam: Elsevier, 2008, pp. 141–169.
[11] K. Janhunen and M.D. Cole. “Development of a predictive model for batch membership of street samples of heroin”. In: Forensic Science International 102.1 (), pp. 1–11.
[12] Patrick Mair, Kurt Hornik, and Jan de Leeuw. “Isotone optimization in R: pool-adjacent-violators algorithm (PAVA) and active set methods”. In: Journal of statistical software 32.5 (2009), pp. 1–24.
[13] Anders Nordgaard, Ricky Ansell, Weine Drotz, and Lars Jaeger. “Scale of conclusions for the value of evidence”. In: Law Probability and Risk 11 (Mar. 2012), pp. 1–24.
[14] John D. Power, John O’Brien, Brian Talbot, Michael Barry, and Pierce Kavanagh. “The identification of an impurity product, 4,6-dimethyl-3,5-diphenylpyridin-2-one in an amphetamine importation seizure, a potential route specific by-product for am-phetamine synthesized by the APAAN to P2P, Leuckart route”. In: Forensic Science In-ternational 241 (), e13–e19.
[15] D. Ramos-Castro and J. González-Rodríguez. “Cross-entropy analysis of the informa-tion in forensic speaker recogniinforma-tion”. In: Odyssey. 2008.
[16] Tania M. G. Salerno, Paola Donato, Giampietro Frison, Luca Zamengo, and Luigi Mon-dello. “Gas Chromatography-Fourier Transform Infrared Spectroscopy for Unambigu-ous Determination of Illicit Drugs: A Proof of Concept”. eng. In: Frontiers in chemistry 8 (July 2020). PMC7396574[pmcid], pp. 624–624.
[17] Scale of conclusions. https://nfc.polisen.se/tjanster/utlatandeskala/. [Online; accessed -23-April-2021]. 2019.
[18] Liang Tang, Silong Peng, Yiming Bi, Peng Shan, and Xiyuan Hu. “A New Method Com-bining LDA and PLS for Dimension Reduction”. In: PLOS ONE 9.5 (May 2014), pp. 1– 10.DOI: 10.1371/journal.pone.0096944.
[19] Liang Tang, Silong Peng, Yiming Bi, Peng Shan, and Xiyuan Hu. “A new method com-bining LDA and PLS for dimension reduction”. In: PloS one 9 (May 2014), e96944. [20] R. Tibshirani, Trevor Hastie, B. Narasimhan, and Gilbert Chu. “Diagnosis of multiple
cancer types by shrunken centroids of gene expression”. In: Proceedings of National Aca-demic Science (PNAS) 99 (Jan. 2002), pp. 6567–6572.
[21] Yukio Tominaga. “Comparative study of class data analysis with PCA-LDA, SIMCA, PLS, ANNs, and k-NN”. In: Chemometrics and Intelligent Laboratory Systems 49.1 (1999), pp. 105–115.
[22] Daan Vangerven. “Development of a semi-automatic chromatogram comparison tool and its application to cannabis resin profiling”. In: (2020).
