# Chapter 7 Trace glass evidence: chemical composition

#### Soyoung Park, Sam Tyner

## 7.1 Introduction

It is easy to imagine a crime scene with glass fragments: a burglar may have broken a glass door, a glass bottle could have been used in an assault, or a domestic disturbance may involve throwing something through a window. During the commission of a crime, there are many ways that glass can break and be transferred from the scene. The study of glass fragments is important to forensic science because the glass broken at the scene can transfer to the perpetrator’s shoes, clothing, or even their hair (Curran, Hicks, and Buckleton 2000).

Crime scene investigators collect fragments of glass at the scene as a part of the evidence collection process, and the fragments are sent to the forensic science lab for processing. Similarly, evidence such as clothing and shoes are collected from a suspect, and if glass is found, the fragments are sent to the lab and compared to the fragments found at the scene. The question that the analyst usually tries to answer is, “Did these glass fragments come from the same source?” This is a source level question, meaning that the comparison of the fragments will only tell the investigators whether or not the fragments from the suspect and the fragments from the scene have the same origin. As discussed in Section 1.3, the forensic analysis will not inform investigators how the suspect came into contact (activity level) with the glass or if the suspect was the perpetrator of the crime (offense level) (Roger Cook et al. 1998).

### 7.1.1 Problems of interest

There are two key problems of interest in glass fragments comparison, but before defining them, we need to define the different glass involved in the investigation of a crime. Glass fragments found on the suspect, for example in their hair, shoes, or clothes, are questioned fragments, which we denote by \(Q\). Glass fragments found at the crime scene, for example in front of a broken window or taken from the broken window, are known fragments, which we denote \(K\). This brings us to a specific source question: Did the questioned fragments \(Q\) found on the suspect come from the same source of glass as the known fragments \(K\), which we know belong to a specific piece of glass at the scene? The goal is now to quantify the similarity between \(Q\) and \(K\). There are lots of ways measure similarity between two glass fragments, but the metric should be defined according to available databases of glass fragment measurements for which ground truth is known. For example, if we have elemental compositions measured in parts per million (ppm) as numerical values, the similarity can be quantified by the difference of the chemical compositions of \(Q\) and \(K\).

### 7.1.2 Current practice

There are many types of glass measurements such as color, thickness, refractive index (RI) and chemical concentrations. In this chapter, we will focus on float glass that is most frequently used in windows, doors and automobiles. For discussions of the other measurements see e.g. Curran, Hicks, and Buckleton (2000). The elemental concentrations of float glass that we use here were obtained through inductively coupled mass spectrometry with a laser add-on (LA-ICP-MS). In the current practice, there are two analysis guides from ASTM International, (ASTM-E2330-12 2012) and (ASTM-E2927-16 2016). To determine the source of glass fragments according to these two guides, intervals around the mean concentrations are computed for each element, and if all elements’ intervals overlap, then the conclusion is that the fragments come from the same source. For more detail on these methods, see ASTM-E2330-12 (2012) and ASTM-E2927-16 (2016).

### 7.1.3 Comparing glass fragments

In order to determine if two glass fragments come from the same source, a forensic analyst considers many properties of the glass, including color, fluorescence, thickness, surface features, curvature, and chemical composition. All methods for examining these properties, except for methods of chemical composition analysis, are non-destructive. If the fragments are large, exclusion are easy to reach if the glass are of different colors because of the wide variety of glass colors possible in manufacturing. Typically, however, glass fragments are quite small and color determination is very difficult. Similarly, thickness of glass is dictated by the manufacturing process, which aims for uniform thickness, so if two glass fragments differ in thickness by more than 0.25mm, an exclusion is made (Bottrell 2009). For glass fragments of the same color and thickness, microscopic techniques for determining light absorption (fluorescence), curvature, surface features (such as coatings), are used before the destructive chemical composition analysis.

### 7.1.4 Goal of this chapter

In this chapter, we construct a new rule for making glass source conclusions using the random forest algorithm to classify the source of glass fragments (Park and Carriquiry 2019a).

## 7.2 Data

### 7.2.1 Chemical composition of glass

The process for determining the chemical composition of a glass fragment is given in great detail in ASTM-E2330-12 (2012) and ASTM-E2927-16 (2016). This destructive method determines elemental composition with Inductively Coupled Plasma Mass Spectrometry (ICP-MS). Up to 40 elements can be detected in a glass fragment using this method. In Weis et al. (2011), only 18 elements are used: calcium (Ca), sodium (NA) and magnesium (Mg) are the major elements, followed by aluminum (Al), potassium (K) and iron (Fe) as minor elements, and lithium (Li), titanium (Ti), manganese (Mn), rubidium (Rb), strontium (Sr), zirconium (Zr), barium (Ba), lanthanum (La), cerium (Ce), neodymium (Nd), hafnium (Hf), and lead (Pb) as the trace elements. The methods of Weis et al. (2011) use standard deviations (\(\sigma\)) of repeated measurements of the same fragment to create intervals around the measurements. Intervals of width \(2\sigma, 4\sigma, 6\sigma, 8\sigma, 10\sigma, 12\sigma, 16\sigma, 20\sigma, 30\sigma,\) and \(40\sigma\) are considered for overlap.

### 7.2.2 Data source

Dr. Alicia Carriquiry of Iowa State University commissioned the collection of a large database of chemical compositions of float glass samples. The details of this database are explained in Park and Carriquiry (2019a). The full database is available here. The database includes 31 panes of float glass manufactured by Company A and 17 panes manufactured by Company B, both located in the United States. The Company A panes are labeled AA, AB, … , AAR, and the Company B panes are labeled BA, BB, … ,BR. The panes from Company A were produced during a three week period (January 3-24, 2017) and the panes from Company B were produced during a two week period (December 5-16, 2016).

To understand variability within a ribbon of glass, two glass panes were collected on almost all days in each company, one from the left side and one from the right side of the ribbon. Twenty four fragments were randomly sampled from each glass pane. Five replicate measurements were obtained for 21 of the 24 fragments in each pane; for the remaining three fragments in each pane, we obtained 20 replicate measurements. Therefore, each pane has 165 measurements for 18 elements. For example, see the heuristic in Figure 7.1. In some panes, there may be a fragment with fewer than five replicate measurements. The unit for all measurements is parts per million (ppm).

    library(tidyverse)
    pane <- expand.grid(x = 1:16, y = 1:4)
    n <- nrow(pane)
    pane$id <- 1:n
    frags <- sample(n, 24)
    rep20 <- sample(frags, 3)
    rep5 <- frags[!(frags %in% rep20)]
    pane_sample <- data.frame(frag = c(rep(rep5, each = 5), rep(rep20, each = 20)),
        rep = c(rep(1:5, 21), rep(1:20, 3)))
    sample_data <- left_join(pane, pane_sample, by = c(id = "frag"))
    ggplot(data = sample_data, aes(x = x, y = y)) + geom_tile(fill = "white", color = "black") +
        geom_jitter(aes(color = as.factor(rep)), alpha = 0.8, size = 0.5) + scale_color_manual(values = c(rep("black",
        20))) + theme_void() + theme(legend.position = "none")

Figure 7.1: An example of how the glass fragments were sampled, if the 64 squares are imagined to be randomly broken fragments within a pane.

### 7.2.3 Data structure

Next, we look at the glass data.

    glass <- read.csv("dat/glass_raw_all.csv")
    head(glass)

    ##   pane fragment Rep mfr element       ppm
    ## 1   AA        1   1   A      Al  2678.000
    ## 2   AA        1   1   A      Ba    10.800
    ## 3   AA        1   1   A      Ca 63140.000
    ## 4   AA        1   1   A      Ce     9.520
    ## 5   AA        1   1   A      Fe   667.000
    ## 6   AA        1   1   A      Hf     1.148

The elements have very different scales, as some (e.g. Ca) are major elements, some (e.g. Al) are minor elements, and others (e.g Hf) are trace elements. Thus, we take the natural log transformation of all measurements to put them on a similar scale.

    # need to make sure the panes are shown in order of mfr date
    pane_order <- c(paste0("A", LETTERS[c(1:13, 15, 22:25)]), paste0("AA", LETTERS[c(1:4,
        6, 8:13, 17:18)]), names(table(glass$pane))[c(32:48)])
    glass$pane <- ordered(glass$pane, levels = pane_order)

    glass_log <- glass %>% mutate(log_ppm = log(ppm))
    glass_log %>% select(-ppm) %>% spread(element, log_ppm) %>% select(mfr, pane,
        fragment, Rep, Li, Na, Mg, Al, K, Ca) %>% head()

    ##   mfr pane fragment Rep        Li       Na       Mg       Al        K
    ## 1   A   AA        1   1 0.8329091 11.54025 10.04715 7.892826 7.136483
    ## 2   A   AA        1   2 0.5423243 11.53772 10.03671 7.883069 7.126248
    ## 3   A   AA        1   3 0.7227060 11.53126 10.04499 7.909122 7.127694
    ## 4   A   AA        1   4 0.7975072 11.54219 10.05449 7.912423 7.144407
    ## 5   A   AA        1   5 0.7227060 11.52968 10.02260 7.871311 7.103322
    ## 6   A   AA        2   1 0.5596158 11.52663 10.04107 7.898411 7.118826
    ##         Ca
    ## 1 11.05311
    ## 2 11.04930
    ## 3 11.04580
    ## 4 11.06147
    ## 5 11.01023
    ## 6 11.03747

    cols <- csafethemes:::csafe_cols_secondary[c(3, 12)]
    glass_log %>% filter(element %in% c("Li", "Na", "Mg", "Al", "K", "Ca")) %>%
        ggplot() + geom_density(aes(x = log_ppm, fill = mfr), alpha = 0.7) + scale_fill_manual(name = "Manufacturer",
        values = cols) + facet_wrap(~element, scales = "free", nrow = 2) + labs(x = "Log concentration (ppm)",
        y = "Density") + theme(legend.position = "top")

Figure 7.2: Density estimation of selected chemical compositions, colored by manufacturers

Figure 7.2 shows density plots of six chemical compositions of Al, Ca, K, Li, Mg, Na. Na and Ca (major elements) show density curves from two manufacturers are overlapped, while Al, K show clear separation between curves by manufacturers. This implies that glass fragments from different manufacturers will be very easy to distinguish.

    glass_log %>% filter(element %in% c("Na", "Ti", "Zr", "Hf")) %>% ggplot() +
        geom_boxplot(aes(x = pane, y = log_ppm, fill = mfr), alpha = 0.8, outlier.size = 0.5,
            size = 0.1) + scale_fill_manual(name = "Manufacturer", values = cols) +
        facet_wrap(~element, scales = "free", nrow = 2, labeller = label_both) +
        theme_bw() + theme(legend.position = "none") + scale_x_discrete(labels = c("AA",
        rep("", 30), "BA", rep("", 15), "BR")) + labs(x = "Pane (in order of manufacture)",
        y = "Log concentration (ppm)")

Figure 7.3: Box plot of four elements in 48 panes, ordered by date of production, from two manufacturers.

Figure @ref(fig:box plot) shows box plots of the measurements of four elements (Na, Ti, Zr, Hf) in each of the 48 panes in the database, colored by manufacturer. Boxes are ordered by date of production within manufacturer. We can confirm the between- and within- pane variability. Interestingly, element values of Zr and Hf in manufacturer A are both decreasing in time, which is evidence that the element measurements are highly correlated. To account for this relationship, we use methods that do not require independence of measurements.

    col_data_AA <- glass_log_long[, -6] %>% filter(pane == "AA") %>% spread(element,
        log_ppm) %>% group_by(fragment) %>% summarise_if(is.numeric, mean, na.rm = TRUE)

    col_data_BA <- glass_log_long[, -6] %>% filter(pane == "BA") %>% spread(element,
        log_ppm) %>% group_by(fragment) %>% summarise_if(is.numeric, mean, na.rm = TRUE)

    P1 <- ggcorr(col_data_AA[, 3:20], geom = "blank", label = TRUE, hjust = 0.75) +
        geom_point(size = 10, aes(color = coefficient > 0, alpha = abs(coefficient) >
            0.5)) + scale_alpha_manual(values = c(`TRUE` = 0.25, `FALSE` = 0)) +
        guides(color = FALSE, alpha = FALSE)

    P2 <- ggcorr(col_data_BA[, 3:20], geom = "blank", label = TRUE, hjust = 0.75) +
        geom_point(size = 10, aes(color = coefficient > 0, alpha = abs(coefficient) >
            0.5)) + scale_alpha_manual(values = c(`TRUE` = 0.25, `FALSE` = 0)) +
        guides(color = FALSE, alpha = FALSE)

    P1 + P2

## 7.3 R Packages

We propose a machine learning method to quantify the similarity between two glass fragments \(Q\) and \(K\). The goal here is to construct a classifier that predicts, with low error, whether or not \(Q\) and \(K\) have the same source. Using the glass database, we construct many pairwise comparisons between two glass fragments with known sources, and we will record their source as either same source or different source.

To construct the set of comparisons, we

1. Take the natural log of all measurements of all glass fragments. (ppm to \(\log\)(ppm))
2. Select one pane of glass in the database to be the “questioned” source. Sample one fragment from this pane, and call it \(Q\).
3. Select one pane of glass in the database to be the “known” source. Sample one fragment from this pane, and call it \(K\).
4. Construct the response variable: if the panes in 2 and 3 are the same, the response variable is KM for known mates. Otherwise, the response variable is KNM for known non-mates.
5. Construct the features: take the mean of the repeated measurements of \(Q\) and \(K\) in each element, and take the absolute value of the difference between the mean of \(Q\) and the mean of \(K\).
6. Repeat 2-5 until we have a data set suitable for training a classifier.

We use the R package `caret` to train a random forest classifier to determine whether two glass fragments (\(Q\) and \(K\)) have the same source or have different sources (Jed Wing et al. 2018).

To begin, the package can be installed from CRAN:

    install.packages("caret")

    library(caret)
    # other package used for plotting
    library(GGally)
    library(patchwork)

The R package `caret` (Classification And REgression Training) is used for applied predictive modeling. The `caret` package allows us to run 238 different models as well as adjusting data splitting, preprocessing, feature selection, tuning parameter, and variable importance estimation. We use it here for fitting a random forest model to our data using cross-validation and down-sampling.

    diff_Q_K_data <- readRDS("dat/rf_data_kfrags_1z.rds")

Table 7.1: Differences of log values of concentrations (Li, Na, Mg, Al, K, Ca) from pairs of known mates (KM) and known non-mates (KNM)
Class | Li | Na | Mg | Al | K | Ca
--- | --- | --- | --- | --- | --- | ---
KM | 0.0078 | 0.0043 | 0.0089 | 0.0270 | 0.0103 | 0.0149
KM | 0.0845 | 0.0065 | 0.0042 | 0.0049 | 0.0126 | 0.0032
KM | 0.0826 | 0.0170 | 0.0069 | 0.0137 | 0.0259 | 0.0014
KNM | 0.5362 | 0.0581 | 0.0267 | 0.0741 | 0.0834 | 0.0619
KNM | 0.2437 | 0.0283 | 0.0092 | 0.0181 | 0.0377 | 0.0386
KNM | 0.2319 | 0.0316 | 0.0283 | 0.0001 | 0.0437 | 0.0518

Table 7.1 shows the example of pairwise differences among glass measurements. If we take the difference of two fragments from the same pane, then KM is assigned to the response variable of `Class`. If we take the difference of two fragments from two different panes, then KNM is assigned to variable `Class`. Each row has 18 differences and one variable `Class` indicating the source of two glass fragments. By taking pairwise differences, there are many more KNM pairs than KM pairs: we can construct 67,680 KNM pairs but only 1,440 KM pairs from the glass database.

    diff_Q_K_data %>% gather(element, diff, Li:Pb) %>% filter(element %in% c("Zr",
        "Li", "Hf", "Ca")) %>% ggplot() + geom_density(aes(x = diff, fill = class),
        alpha = 0.7) + scale_fill_manual(name = "Class", values = cols) + facet_wrap(~element,
        scales = "free", nrow = 2) + labs(x = "Difference between Q and K (log(ppm))") +
        theme(legend.position = "top")

Figure 7.4: Histogram of differences from four chemical elements among KM and KNM.

Figure 7.4 shows the distribution of differences for KM and KNM pairs in the Ca, Hf, Li, and Zr. Across all elements, distribution of differences from KNM pairs are more dispersed than differences of KM. The differences of KM have high density near zero, while the KNM measurements are shifted to the right and have a long tail. These differences for all 18 elements will be the features of the random forest.

Finally, we construct the random forest classifier. Since the data have 1,440 KM and 67,680 KNM observations, the response variable is imbalanced. With this large imbalance, the algorithm can simply predict KNM and have low error rate without learning anything about the properties of the KM class. You can find more ways to consider this imbalance problem for fitting the random forest in this glass fragment source prediction in (Park and Carriquiry 2019a). In this chapter, we will down-sample the KNM observations to equal the number of KM observations. That way, we will have 1,440 KM and 1,440 KNM comparisons. Then, we sample 70% of them to be the training set and the remaining 30% are the testing set.

    diff_Q_K_data$class <- as.factor(diff_Q_K_data$class)
    # Down sample the majority class to the same number of minority class
    set.seed(123789)
    down_diff_Q_K_data <- downSample(diff_Q_K_data[, c(1:18)], diff_Q_K_data$class)
    down_diff_Q_K_data <- down_diff_Q_K_data %>% mutate(id = row_number())
    names(down_diff_Q_K_data)[19] <- "class"
    table(down_diff_Q_K_data$class)

    ##
    ##   KM  KNM
    ## 1440 1440

    # Create training set
    train_data <- down_diff_Q_K_data %>% sample_frac(0.7)
    # Create test set
    test_data <- anti_join(down_diff_Q_K_data, train_data, by = "id")
    train_data <- train_data[, -20]  #exclude id
    test_data <- test_data[, -20]  #exclude id
    # dim(train_data) # 2016 19 dim(test_data) # 864 19

After down-sampling and setting aside 30% of the data for testing, we can train the classifier. Below is the R code to fit the random forest, using `caret` (Jed Wing et al. 2018). The tuning parameter for the random forest algorithm is `mtry`, the number of variables available for splitting at each tree node. We use five values (1, 2, 3, 4, 5) to tune `mtry` and pick the optimal one. For the classification, the default `mtry` is the square root of the number of predictor variables (it is \(\sqrt{18} \approx 4.2\), in our study). To pick the optimal `mtry` parameter, the area under the receiver operating characteristic (ROC) curve is calculated. The optimal value will result the highest area under the ROC curve (AUC). The data are also automatically centered and scaled for each predictor. We use 10-fold cross-valildation to evaluate the random forest algorithm and repeat entire process three times.

    ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3, savePredictions = "final",
        summaryFunction = twoClassSummary, classProbs = TRUE)

    # mtry is recommended to use sqrt(# of variables)=sqrt(18) so 1:5 are tried
    # to find the optimal one
    RF_classifier <- train(class ~ ., train_data, method = "rf", tuneGrid = expand.grid(.mtry = c(1:5)),
        metric = "ROC", preProc = c("center", "scale"), trControl = ctrl)

    RF_classifier <- readRDS("dat/RF_classifier.RDS")
    RF_classifier

    ## Random Forest
    ##
    ## 2016 samples
    ##   18 predictor
    ##    2 classes: 'KM', 'KNM'
    ##
    ## Pre-processing: centered (18), scaled (18)
    ## Resampling: Cross-Validated (10 fold, repeated 3 times)
    ## Summary of sample sizes: 1815, 1814, 1814, 1814, 1815, 1814, ...
    ## Resampling results across tuning parameters:
    ##
    ##   mtry  ROC        Sens       Spec
    ##   1     0.9637918  0.9826797  0.8323266
    ##   2     0.9641196  0.9807190  0.8467239
    ##   3     0.9625203  0.9758170  0.8457239
    ##   4     0.9624018  0.9738562  0.8467172
    ##   5     0.9614865  0.9728758  0.8473838
    ##
    ## ROC was used to select the optimal model using the largest value.
    ## The final value used for the model was mtry = 2.

From the RF fitting result, it showed the AUC (ROC), Sens (Sensitivity) and Spec (Specificity) values in each mtry values. R-package caret selects the best mtry value (it is 2, in our study) to give the highest AUC (ROC) value 0.964, from 10-fold cross-validation and 3 repeated process.

The final performance of the random forest algorithm with the optimal mtry of 2 with 500 of number of trees. In the training set, the false negative rate is 0.02 and the false positive rate is 0.153 of false positive rate, which is the optimized classification result. We see the OOB estimate of error rate of 0.085, which is the average error rate from the fold left out of the training set during cross-validation.

    RF_classifier$finalModel$confusion

    ##       KM KNM class.error
    ## KM  1000  20  0.01960784
    ## KNM  152 844  0.15261044

    imp <- varImp(RF_classifier)$importance
    imp <- as.data.frame(imp)
    imp$varnames <- rownames(imp)  # row names to column
    rownames(imp) <- NULL
    imp <- imp %>% arrange(-Overall)
    imp$varnames <- as.character(imp$varnames)
    elements <- read_csv("dat/elements.csv")

    imp %>% left_join(elements[, c(2, 4)], by = c(varnames = "symb")) %>% ggplot(aes(x = reorder(varnames,
        Overall), y = Overall, fill = classification)) + geom_bar(stat = "identity",
        color = "grey40") + scale_fill_brewer(name = "Element", palette = "Blues",
        direction = -1) + labs(y = "Overall importance", x = "Variable") + scale_y_continuous(position = "right") +
        coord_flip() + theme_bw()

Figure 7.5: Variable importance from the RF classifier, colored by element types (major, minor or trace).

By fitting the random forest algorithm, we can also get a measure of variable importance. This metric ranks which elements out of 18 are most important to correctly predicting the source of the glass fragments. Figure 7.5 shows that elements K, Ce, Zr, Rb, and Hf are five very important variables and Pb, Sr, Na, Ca, and Mg are five minimally important variables. All major elements are not important: this is not surprising because all glass is, broadly speaking, very similar chemically. Conversely, most of the trace elements are ranked in the top half.

## 7.4 Drawing Conclusions

The result of the random forest classifier is a probability of the observation belonging to the two classes, KM and KNM, and a predicted class. For an observation, if the class probability for KM is greater than 0.5, then the class prediction is M (mated), otherwise it is NM (non-mated). We use the random forest classifier to predict on the test set. Recall that we know the ground truth, KM and KNM for the test data.

    # Get prediction from RF classifier in test data
    pred_class <- predict(RF_classifier, test_data)
    table_pred_class <- confusionMatrix(pred_class, test_data$class)$table
    table2 <- data.frame(table_pred_class) %>% spread(Reference, Freq)
    table2$Prediction <- c("M", "NM")
    table2 %>% knitr::kable(format = "html", caption = "Classification result of test set",
        longtable = FALSE, row.names = FALSE)

Table 7.2: Classification result of test set
Prediction | KM | KNM
--- | --- | ---
M | 409 | 69
NM | 11 | 375

Table 7.2 shows the classification results on the testing set from the random forest. There are 420 KM and 444 KNM comparisons in the test set. For the KM cases, the RF classifier correctly predicts the source as M in 97.4% of cases and wrongly predicts the source as NM in 2.6% of cases. This is the false negative rate, FNR. For the KNM cases, the RF correctly classifies them as NM 84.5% of the time, while 15.5% are incorrectly classified as M. We know that the method resulted in higher FPR (15.5) because the data we have are from many close panes made by the same manufacturer that are difficult to distinguish from one another.

    # class probability for the same pane is used as similarity score for RF
    prob_prediction <- predict(RF_classifier, test_data, type = "prob")[, "KM"]
    test_data$prob <- prob_prediction

    test_data %>% ggplot() + geom_density(aes(x = prob, fill = class), alpha = 0.7) +
        scale_fill_manual(name = "Class", values = cols) + labs(x = "Empirical probability for KM from RF classifier",
        y = "Density") + theme(legend.position = "top")

Figure 7.6: Scores from the RF classifier for the test set. Ground truth (either KM or KNM) is known. Any observations above 0.5 are declared same source, while all others are declared different source. Notice the long tail in the KNM cases.

Figure 7.6 shows the distribution of the RF scores, colored by true classes in test data. The score is the empirical class probability that an observation belongs to the KM class. The modes of the the two classes are well separated, while there is still some overlap between the classes’ density curve. The tail of the distribution of scores for KNMs is more skewed than that of the KMs, which helps explain the higher false positive rate. Table 7.2 is the predicted classification result using the cut-off of 0.5. If the RF score is larger than 0.5, then we will predict the pair of glass fragments have the same source (M). If not, then we declare they have different source (NM).

## 7.5 Case Study

Here we introduce a new set of five pairs of comparisons and use the random forest classifier we trained in Section @ref(glass_rpkgs) to determine if the two samples being compared are from the same source of glass or from different sources.

Suppose you are given the following data on eight glass samples, where each sample has been measured 5 times:

Your supervisor wants you to make five comparisons:

compare_id | sample_id_1 | sample_id_2
--- | --- | ---
1 | 1 | 2
2 | 1 | 3
3 | 1 | 5
4 | 4 | 6
5 | 7 | 8

First, you need to take the log of the measurements, then get the mean for each sample and each element.

    # take the log, then the mean of each for comparison
    new_samples <- mutate(new_samples, logppm = log(ppm)) %>%
      select(-ppm) %>%
      group_by(sample_id, element) %>%
      summarise(mean_logppm = mean(logppm))

Next, make each of the five comparisons. We write a function, `make_compare()` to do so, then we use `purrr::map2()` to perform the comparison for each row in `compare_id`. The resulting compared data is below.

    make_compare <- function(id1, id2) {
        dat1 <- new_samples %>% filter(sample_id == id1) %>% ungroup() %>% select(element,
            mean_logppm)
        dat2 <- new_samples %>% filter(sample_id == id2) %>% ungroup() %>% ungroup() %>%
            select(element, mean_logppm)
        tibble(element = dat1$element, diff = abs(dat1$mean_logppm - dat2$mean_logppm)) %>%
            spread(element, diff)
    }

    compared <- comparisons %>% mutate(compare = map2(sample_id_1, sample_id_2,
        make_compare)) %>% unnest()

    DT::datatable(compared) %>% DT::formatRound(4:21, 3)

Finally, we take the compared data and predict using the random forest object:

    newpred <- predict(RF_classifier, newdata = compared, type = "prob")

    # only take the id data from compare, and the M prediction from newpred
    bind_cols(compared[, 1:3], score = newpred[, 1]) %>% mutate(prediction = ifelse(score >
        0.5, "M", "NM")) %>% knitr::kable(format = "html", caption = "Predicted scores of new comparison cases from the random forest classifier",
        longtable = FALSE, row.names = FALSE, digits = 3, col.names = c("Comparison ID",
            "Sample ID 1", "Sample ID 2", "Score", "Decision"))

Table 7.3: Predicted scores of new comparison cases from the random forest classifier
Comparison ID | Sample ID 1 | Sample ID 2 | Score | Decision
--- | --- | --- | --- | ---
1 | 1 | 2 | 0.812 | M
2 | 1 | 3 | 0.238 | NM
3 | 1 | 5 | 0.002 | NM
4 | 4 | 6 | 0.430 | NM
5 | 7 | 8 | 0.618 | M

Table 7.3 shows the RF score which is defined as the empirical class probability that the samples are from the same pane (M) by the RF algorithm6. Based on these scores, it appears that samples 1 and 2 are very likely from the same source, while samples 1 and 5 are very likely from different sources. Using a threshold of 0.5, sample 4 and sample 6 are predicted to be from different sources. However, this score is so close to the threshold of 0.5 that we may want to say this is inconclusive.7 Samples 1 and 3 are probably from different sources, though we are less confident in that determination than we are in the decision that samples 1 and 5 are from different sources. Finally, samples 7 and 8 are probably from the same source, but we are less certain of this than we are that samples 1 and 2 are from the same source.

More research into the appropriateness of this random forest method for making glass source conclusions is needed. Important points to consider are:

- What database should be used to train the algorithm?
- Could the random forest method over fit, reducing the generalizability of this method?

Nevertheless, the RF classifier makes good use of the glass fragment data, with its high dimension and small number of repeated measurements. More details on the RF classifier and additional discussion can be found in Park and Carriquiry (2019a).

### References

Curran, James Michael, Tacha Natalie Hicks, and John S. Buckleton. 2000. Forensic Interpretation of Glass Evidence. CRC Press.

Cook, Roger, Ian W Evett, Graham Jackson, PJ Jones, and JA Lambert. 1998. “A Hierarchy of Propositions: Deciding Which Level to Address in Casework.” Science & Justice 38 (4): 231–39.

ASTM-E2330-12. 2012. “Standard Test Method for Determination of Concentrations of Elements in Glass Samples Using Inductively Coupled Plasma Mass Spectrometry (ICP-MS) for Forensic Comparisons, Retrieved from Https://Doi.org/10.1520/E2330-12.” ASTM International.

ASTM-E2927-16. 2016. “Standard Test Method for Determination of Trace Elements in Soda-Lime Glass Samples Using Laser Ablation Inductively Coupled Plasma Mass Spectrometry for Forensic Comparisons, Retrieved from Https://Doi.org/10.1520/E2927-16.” ASTM International.

Bottrell, Maureen C. 2009. “Forensic Glass Comparison: Background Information Used in Data Interpretation. (Report).” Forensic Science Communications 11 (2). https://archives.fbi.gov/archives/about-us/lab/forensic-science-communications/fsc/april2009/review.

Park, Soyoung, and Alicia Carriquiry. 2019a. “Learning Algorithms to Evaluate Forensic Glass Evidence.” Annals of Applied Statistics. https://www.e-publications.org/ims/submission/AOAS/user/submissionFile/36724?confirm=02d8a608.

Weis, Peter, Marc Dücking, Peter Watzke, Sonja Menges, and Stefan Becker. 2011. “Establishing a match criterion in forensic comparison analysis of float glass using laser ablation inductively coupled plasma mass spectrometry.” Journal of Analytical Atomic Spectrometry 26 (6): 1273–84. https://doi.org/10.1039/c0ja00168f.

Jed Wing, Max Kuhn. Contributions from, Steve Weston, Andre Williams, Chris Keefer, Allan Engelhardt, Tony Cooper, Zachary Mayer, et al. 2018. Caret: Classification and Regression Training. https://CRAN.R-project.org/package=caret.

6. In fact, we do know the truth here: for comparisons (1, 2, 3, 4, 5), the truth is (M, NM, NM, M, M).↩

7. See this discussion of “equivocal zones:” http://appliedpredictivemodeling.com/blog/2013/8/15/equivocal-zones↩