Enhancing the classification of isolated theropod teeth using machine learning: a comparative study

Data processing and preparation for the modelling phase

The pre-processing phase is divided into three different stages: the first corresponds to data cleaning, which includes selection of variables and taxa, the second corresponds to the oversampling of the data, and the third corresponds to the standardization of the data.

The theropod teeth measurements were obtained from Hendrickx et al. (2023) (Fig. 1), which is an updated dataset with morphometric variables on a large set of theropod teeth compiled from previous works (e.g., Smith, Vann & Dodson, 2005; Hendrickx, Mateus & Araújo, 2015; Hendrickx, Tschopp & Ezcurra, 2020). The variables initially evaluated were: (1) Crown Base Length (CBL); (2) Crown Base Width (CBW); (3) Crown Height (CH); (4) Apical Length (AL); (5) Crown Base Ratio (CBR); (6) Crown Height Ratio (CHR); (7) Mid-crown Length (MCL); (8) Mid-crown Width (MCW); (9) Mid-crown Ratio (MCR); (10) Mesiobasal Denticle Extension (MDE); (11) Mesial Serrated Carina Length (MSL); (12) Crown Angle (CA); (13) Mesioapical Denticle Density (MA); (14) Mesiocentral Denticle Density (MC); (15) Mesiobasal Denticle Density (MB); (16) Distoapical Denticle Density (DA); (17) Distocentral Denticle Density (DC); (18) Distobasal Denticle Density (DB), Distal Denticle Length (DDL); (19) Average Mesial Denticle Density (MAVG); (20) Average Distal Denticle Density (DAVG); and (21) Denticle Size Density Index (DSDI). The isolated theropod teeth were found mainly in the United States (U.S.- 435 observations), Canada (300 observations), Argentina (99 observations) and China (89 observations) (Table 1). The ontogenetic state of the individuals is mostly uncertain (Table 2), with most specimens belonging to adult and immature individuals in different ontogenetic stages. The geological epoch of the deposits where the teeth were collected was variate, with the most common epochs belonging to Late Cretaceous (724 observations), “Middle” Cretaceous (238 observations), and Late Jurassic (205 observations) (Table 3). For more information about the methodology for the measurement of the variables and information on the phylogenetic framework of the isolated theropod teeth used, see Hendrickx, Mateus & Araújo (2015).

Data pre-processing began by removing variables that contained more than 15% missing values from the total number of isolated teeth, this resulted in the selection of the following variables: DC, DDL, CBW, CBL, CBR, CHR and CH (Fig. 2). The inclusion of both raw variables (CBW, CBL, CH) and their ratio-derived counterparts (CBR and CHR) in the analysis was chosen since most machine learning methods, such as gradient boosting trees and random forests are less sensitive to multicollinearity as they prioritize variable splits rather than linear dependencies. By keeping both types of variables, we increase the possibilities of the models identifying subtle and crucial patterns in the dataset that might otherwise be missed if ratios or raw variables were excluded. After data cleaning, two datasets were created. The first for classification at the genus level, and the second for classification at a higher taxonomic level. The second dataset allowed for the inclusion of groups with fewer data cases that were excluded from the genus-level dataset to avoid class imbalance.

To maintain dataset quality, taxa within each dataset that contained observations equal to or fewer than the number of selected variables were removed, following the methodology of Wills, Underwood & Barrett (2020). After this data processing, the same method of removal was applied to the previously excluded taxonomic groups, but removing the information for the tooth position, which was done to keep as much information as possible. After these steps, the higher taxonomic-level dataset contained 1,119 isolated theropod teeth belonging to twenty-three combinations of taxa and tooth position, including Abelisauridae lateral, Abelisauridae mesial, Allosauridae lateral, Allosauridae mesial, early-branching Theropoda, Carcharodontosauridae lateral, Carcharodontosauridae mesial, Dromaeosauridae lateral, Dromaeosauridae mesial, Metriacanthosauridae lateral, Neovenatoridae lateral, Noasauridae lateral, Non-abelisauroid Ceratosauria lateral, Non-abelisauroid Ceratosauria mesial, Non-averostran Neotheropoda lateral, Non-spinosaurid Megalosauroidea lateral, Non-tyrannosaurid Tyrannosauroidea lateral, Non-tyrannosaurid Tyrannosauroidea mesial, Spinosauridae, Therizinosauria, Troodontidae lateral, Tyrannosauridae lateral, Tyrannosauridae mesial; and the genus-level dataset contained 925 isolated theropod teeth from thirty-seven combinations of genera and tooth position, including Acrocanthosaurus lateral, Albertosaurus lateral, Alioramus lateral, Allosaurus lateral, Allosaurus mesial, Atrociraptor lateral, Australovenator, Bambiraptor lateral, Baryonyx lateral, Carcharodontosaurus lateral, Ceratosaurus lateral, Ceratosaurus mesial, Coelophysis lateral, Daspletosaurus lateral, Deinonychus lateral, Dilong lateral, Dromaeosaurus lateral, Eoraptor, Fukuiraptor lateral, Gorgosaurus lateral, Majungasaurus lateral, Majungasaurus mesial, Masiakasaurus lateral, Monolophosaurus, Pectinodon lateral, Raptorex lateral, Richardoestesia lateral, Saurornitholestes lateral, Sinraptor, Suchomimus lateral, Torvosaurus lateral, Troodon lateral, Tyrannosaurus lateral, Tyrannosaurus mesial, Velociraptor lateral, Zapsalis lateral.

To obtain more reliable estimate, k-fold cross-validation was used to generate the training and testing datasets, as has been done in other works (Wills, Underwood & Barrett, 2020; Hendrickx et al., 2023). The k-fold cross-validation divides the dataset into k subsets (folds), with each fold used as a testing set while the remaining k-1 folds are used for training. Following the methodology used by Hendrickx et al. (2023), a value of five was chosen for k. Stratified k-fold cross-validation was used instead of standard k-fold cross-validation to ensure that each fold preserves the class distribution within the dataset, preventing loss of information from underrepresented classes. After splitting the data using cross-validation, oversampling and standardization was applied to the training dataset, and the same standardization parameters were subsequently applied to the testing dataset. The different methods of oversampling were applied to avoid possible model bias for more represented taxa. The following four oversampling techniques used in previous studies (e.g., Nguyen, Cooper & Kamei, 2011; Douzas, Bacao & Last, 2018; Wills, Underwood & Barrett, 2020) were evaluated: (1) Random Oversampling (RandomOverSampler), which repeats values (Fernández et al., 2018) and is a non-heuristic algorithm; (2) SMOTE, which generates synthetic samples for the minority class by interpolating between existing minority samples (Chawla et al., 2002); (3) Borderline Synthetic Minority Oversampling Technique (BorderlineSMOTE), which is a variant of SMOTE, but focuses on generating synthetic samples near the borderline of the classes (Han, Wang & Mao, 2005); and (4) K-Means Synthetic Minority Oversampling Technique (KMeansSMOTE), which uses K-Means clustering to generate synthetic samples by clustering the minority class and then applying SMOTE within each cluster (Douzas, Bacao & Last, 2018). We then compared three scaling methods, which were chosen following Ahsan et al. (2021) that analyzes the effect of different standardization techniques on diverse ML model performances. Firstly, we used the standard scaler (StandardScaler), which transforms the data so that the mean is 0 and the standard deviation (SD) is 1 (Sola & Sevilla, 1997). It uses the Z-score normalization, where the standardized value is given by: (1)${\displaystyle z=\frac{\left(x-\mu \right)}{\sigma }}$

where x is the original value, µthe mean value of the variable, and σ is the SD of the variable.

Secondly, we used the robust scaler (RobustScaler), which scales the data according to the inter-quartile range (IQR) (Aggarwal, 2017), where the standardized value is given by: (2)${\displaystyle x^{\prime }=\frac{x-x_{median}}{Q_3-Q_1}}$

where x_median is the median value of the variable, and Q1 and Q3 refer to the first and third quartiles, respectively.

Finally, we used the Quantile Transformer (QuantileTransformer), which transforms the data to follow a normal distribution, with values between 0 and 1 (Pedregosa et al., 2011), where the standardized value is obtained by first getting the value of: (3)${\displaystyle x^{\prime }={\Phi }^{-1}\left(\frac{\text{rank}\left(x\right)}{N+1}\right)}$

where Φ⁻¹ is the inverse of the cumulative distribution function of the standard normal distribution, rank(x) is the rank of the value x in the sorted data, and N is the total number of data points.

Data modeling

Using the training and testing datasets, we evaluated the performance of the following six machine learning multi-class classification models: (1) RF, which was previously used in theropods teeth classification (Wills, Underwood & Barrett, 2020; Wills, Underwood & Barrett, 2023; Hendrickx et al., 2023); (2) gradient boosting (GB), which from our knowledge has not yet been used on theropod studies, but has been applied in different dentistry problems (Carrillo-Perez et al., 2021); (3) K-NN and (4) SVM, which were used to distinguish different types of dinosaur hunters (Millar, 2019); (5) NN, which was previously used to infer dinosaur diets based on dental morphology (Ballell, Benton & Rayfield, 2022); and (6) discriminant analysis, which has been commonly used on theropods teeth classification (e.g., Wills, Underwood & Barrett, 2020; Wills, Underwood & Barrett, 2023; Hendrickx et al., 2023). RF is an ensemble method that uses multiple decision trees (Breiman, 2001). This model uses Gini index to measure the quality of a split. For example, if a target is a classification outcome taking on values 0, 1, …, k − 1, for node m, the proportion of class k observations in node m is given by: (4)${\displaystyle p_{m,k}=\frac{1}{n_m}\sum _{y\in Q_m}I\left(y=k\right)}$

The measure of the impurity is given by ‘Gini’ Yuan, Wu & Zhang (2021) which is obtained by: (5)${\displaystyle H\left(Q_m\right)=\sum _k{p}_{mk}\left(1-p_{mk}\right)}$

where Q_m represents the data node at node m, and n_m is the number of samples. The minimum number of samples required to split an internal node is 2. The minimum number of samples required to be at a leaf node is 1 and the maximum depth of the decision trees was 22. GB is also an ensemble technique that builds trees sequentially, where each tree corrects errors of the previous ones (Friedman, 2001). We used a learning rate of 0.1 and the log-loss function (Murphy, 2012), is given by: (6)${\displaystyle \text{Log Loss}=-\left[y\cdot log\left(p\right)+\left(1-y\right)\cdot log\left(1-p\right)\right]}$

where y is the true label, p is the predicted probability for a given class. The number of boosting stages to perform was 100, the fraction of samples used for fitting the individual base learners was 1, and the function to measure the quality of a split was the mean squared error (MSE) with improvement score by Friedman (2001), which is given by: (7)${\displaystyle \text{Friedman MSE}=\frac{1}{N}\sum _{i=1}^N\left({\left(y_i-{\hat{y}}_i\right)}^2\right)-\lambda \cdot \frac{1}{N}\sum _{i=1}^N{\hat{y}}_i^2}$

where N is the number of samples, y_i is the actual target value for the ith sample, and λ is a regularization parameter that controls the trade-off between variance and bias. K-NN method classifies a sample based on the majority class of its K nearest neighbors (Cover & Hart, 1967). We used K = 5, and all points in each neighborhood are weighted equally, following a uniform distribution. The metric used for the distance computation was the Minkowsk distance (Jain, Duin & Mao, 2000), which is obtained by: (8)${\displaystyle \parallel u-v{\parallel }_2={\left(\sum {\left|u_i-v_i\right|}^p\right)}^{1/p}}$

where u and v are input arrays, p is the order of the norm of the difference |u − v|. NN method uses layers of neurons to learn (Goodfellow, Bengio & Courville, 2016). The models’ hyper-parameters were determined through preliminary runs. The model was trained for 50 epochs of training, the number of training examples used in one forward and backward pass through the network (batch size) was 32. The network architecture comprised three fully connected layers: a first layer containing 128 neurons, followed by a second layer with 64 neurons, and a third layer with 32 neurons, all these layers with a Rectified Linear Unit (ReLU) activation function. The output layer was a Softmax layer, where the output size corresponded to the number of target classes in the dataset. The ReLU activation function (Ying et al., 2019) value is obtained by: (9)${\displaystyle \text{ReLU}\left(x\right)=max\left(0,x\right)}$

The optimizer used was the Adaptive Moment Estimation (Adam) Kingma & Ba (2014) which uses the following initialize moment vectors: m₍₀₎ = 0 and v₍₀₎ = 0 , where m_t and v_(t) are estimates of the first moment (mean) and the second moment (uncentered variance) of the gradients respectively. The update biased first moment estimate is given by: (10)${\displaystyle m_t={\beta }_1\cdot m_{t-1}+\left(1-{\beta }_1\right)\cdot g_t}$

and the update biased second moment estimate is given by: (11)${\displaystyle v_t={\beta }_2\cdot v_{t-1}+\left(1-{\beta }_2\right)\cdot g_t^2}$

where g_t is the gradient at step t, β₁ and β₂ are the decay rates for the moment estimates. To compute the bias-corrected first moment estimate the following formula is used: (12)${\displaystyle {\hat{m}}_t=\frac{m_t}{1-{\beta }_1^t}}$

To compute the bias-corrected second moment estimate the following formula: (13)${\displaystyle {\hat{v}}_t=\frac{v_t}{1-{\beta }_2^t}}$

The parameter update rule is given by: (14)${\displaystyle {\theta }_t={\theta }_{t-1}-\alpha \cdot \frac{{\hat{m}}_t}{\sqrt{{\hat{v}}_t}+ϵ}}$

where θ_t represents the parameters being optimized, ϵ is a small constant to prevent division by zero and α is the learning rate.

The loss function used was the Sparse Categorical Cross-Entropy (Abadi et al., 2016) and its values is given by: (15)${\displaystyle \text{Loss}=-\frac{1}{N}\sum _{i=1}^{N}log\left(p\left(y_i\right)\right)}$

where p(y_i) is the predicted probability of the true class y_i, y_i is the true class label for the ith sample and N is the number of samples. The Softmax Bridle (1989) value is given by: (16)${\displaystyle \sigma \left(z_i\right)=\frac{e^{z_i}}{\sum _{j=1}^K{e}^{z_j}}}$

where z is a vector with the raw prediction scores, K is the number of classes, e^z_i is the exponential of the score for class i, ${\sum }_{j=1}^K{e}^{z_j}$ is the sum of exponentials of all the scores across all classes, and σ(z_i) represents the probability of the input belonging to class i. The model learns after each epoch of training while minimizing the loss value. Quadratic discriminant analysis (QDA) model’s class-specific covariance matrices allowing the formation of quadratic decision boundaries, where the class proportions are inferred from the training data and by fits class conditional densities to the data and using Bayes’ rule (Friedman, 1989), given by: (17)${\displaystyle P\left(C_k\mid x\right)=\frac{P\left(x\mid C_k\right)P\left(C_k\right)}{P\left(x\right)}}$

where P(C_k∣x) is the posterior probability of class k given the feature vector x, P(x∣C_k) is the likelihood of observing x given that it belongs to class k, P(C_k) is the prior probability of class k, P(x) is the marginal likelihood of x. The absolute threshold for a singular value to be considered significant is 1.0e−4. The final classification decision is based on the equation: (18)${\displaystyle \hat{C}\left(x\right)=argmax\left(\frac{1}{2}log\left|{\Sigma }_k\right|-\frac{1}{2}{\left(x-{\mu }_k\right)}^T{\Sigma }_k^{-1}\left(x-{\mu }_k\right)+logP\left(C_k\right)\right)}$

where μ_k is the mean vector for class k, Σ_k is the covariance matrix for class k and x is the feature vector. A regularization parameter (λ) was included with a value of 0.1. This parameter regularizes the per-class covariance estimates by transforming the scaling attribute of a given class (S²) as: (19)${\displaystyle S^2=\left(1-\lambda \right)\Sigma +\lambda I_n}$

where Σ is the estimated covariance matrix of the data, n corresponds to the number of features and I_n to the identity matrix. Finally, SVM method finds the optimal hyperplane to separate different classes in the feature space (Cortes & Vapnik, 1995). The kernel type used in the algorithm was the radial basis function kernel (rbf) (Ding et al., 2021), which is given by: (20)${\displaystyle K\left(x_1,x_2\right)=exp\left(-\gamma \cdot \parallel x_1-x_2{\parallel }^2\right)}$

where γ is the influence of each individual training sample on the decision boundary, x₁ and x₂ are two data points and ∥x₁ − x₂∥² is the Euclidean distance between the two points. The degree of the polynomial kernel function is 3.

A total of 120 model combinations (classifier, oversampling and standardization methods) were trained and evaluated. Since we used two different classification tasks and used two different variable types (log-scaled and non-log-scaled), we had a total of 480 models trained, evaluated and compared.

Model evaluation metrics

All models used a random number generator, and a fixed seed was used for reproducibility. The influence of the variables on the classification results was assessed using permutation feature importance. This technique evaluates the significance of each variable in a predictive model by randomly shuffling the values of a single feature, thereby breaking its relationship with the target. The resulting degradation in the model’s performance indicates how much the model relies on that particular feature.

All models were compared using different metrics, including accuracy (correct classification rate), precision (positive predictive value), and F1-score. These metrics are given by: (21)${\displaystyle \text{Accuracy}=\frac{TP+TN}{TP+TN+FP+FN}}$ (22)${\displaystyle \text{Precision}=\frac{TP}{TP+FP}}$ (23)${\displaystyle \text{Recall}=\frac{TP}{TP+FN}}$ (24)${\displaystyle \text{F1-Score}=\frac{2\times \text{Precision}\times \text{Recall}}{\text{Precision}+\text{Recall}}}$

where TP represents the correctly predicted positive cases, TN represents the correctly predicted negative cases, FP represents the cases incorrectly predicted as positive, and FN represents the cases incorrectly predicted as negative.

Comparison between log-scaled and non-log-scaled variables

Log-scaling is a data transformation technique used to address skewed distributions and reduce the impact of outliers. By compressing the range of values, particularly in data spanning multiple orders of magnitude, this transformation can reveal underlying patterns and potentially improve model performance (Iglewicz & Hoaglin, 1993). Log-scaling has been commonly applied in morphometric analyzes of isolated theropod teeth (e.g., Malafaia et al., 2017; Hendrickx et al., 2019; Wills, Underwood & Barrett, 2020; Wills, Underwood & Barrett, 2023; Hendrickx et al., 2023). To assess the effect of this transformation, a comparison was made between models trained on non-log-scaled and log-scaled variables to evaluate its impact on model performance. This comparison was done by computing the difference between the mean metrics across all the testing folds using the non-log-scaled and log-scaled variables. Log-scaling was done using the natural logarithm.

Implementation information

All data pre-processing, including variables and taxa selection, as well as the correlation analyzes, was performed with the R Software (R Core Team, 2022), version 4.2.1. All standardizations and oversampling methods were applied using Python (Python Software Foundation, 2024), version 3.8.16, using scikit-learn (Pedregosa et al., 2011) and imbalanced-learn (Lemaître, Nogueira & Aridas, 2017) respectively. All metrics values and permutation results were obtained using scikit-learn (Pedregosa et al., 2011). All the models (except the NNs) were created using scikit-learn (Pedregosa et al., 2011), and the NN models were created using TensorFlow (Abadi et al., 2016). The trained models and scalers corresponding to the best models for each classification task can be found in Github. These trained models can be used directly to classify the probability of new isolated theropod teeth belonging to each class.

Application of the methodology in the sample of isolated theropod teeth from Guimarota fossil site

The classification methods were also tested on a sample of unpublished theropod teeth (except MG 27808_D118 that was studied by Zinke (1998) based on morphological comparisons) collected during the last decades of the 20th century in the Kimmeridgian (Upper Jurassic) Guimarota fossil site (Leiria, Portugal). The initial classification here mentioned is based on the labels associated with the specimens, which probably correspond to preliminary identifications made by German researchers that first processed the fossil record from this locality (Table 4 and Fig. 3). These specimens are currently housed by Museu Geológico of Laboratório Nacional de Engenharia e Geologia (LNEG- Lisboa).

Classification at genus level

The original dataset contained information of 1,371 isolated theropod teeth belonging to 93 genera. The dataset after data processing includes 909 isolated theropod teeth from thirty-six combinations of genera and tooth positions. There are two predominant genus combinations: Saurornitholestes lateral (133 observations) and Tyrannosaurus lateral (106 observations) (Fig. 4). After applying oversampling techniques, all training classes were balanced to contain the same number of observations each.

The oversampling influence on the models varied between the models. While GB and SVM obtained lower evaluation metrics when oversampling methods were not applied, the K-NN, NN and QDA classifiers presented higher metric results when not using oversampling methods (see Fig. 5 for F1-score and the ShinyApp for the other metrics results). The model performance varied significantly across oversampling and scaling techniques. RF and GB were the best performing models, with consistent results between different oversampling and standardization techniques. GB benefited from RandomOversampler, K-NN and NN improved their performance by using standardization methods (e.g., StandardScaler and RobustScaler), but they remained less effective when compared to RF and GB. QDA and SVM presented the worst overall metric results and an increase in the metric results when using standardization techniques, namely with StandardScaler and QuantileTransformer (see Fig. 5 for F1-score and the ShinyApp for the other metrics results).

For the non-log-scaled variables, RF and GB were the top-performing classifiers, consistently achieving high metric values across different oversampling and standardization approaches, with average performance metrics ranging from 67.3% to 75.6% on the test set and consistently 100% on the train set. The SVM and QDA models produced the worst results, with metric values ranging from 64.1% to 74.1% on the test and from 74.3% to 80.9% on the train set. Among the tested oversampling methods, BorderlineSMOTE and SMOTE achieved the highest metrics, with testing results ranging from 66.8% to 73.4% across different models and standardization approaches. In contrast, NoOversampling and KMeansSMOTE had the worst performance, with metric values ranging from 66.6% to 73.7%. Regarding standardization techniques, StandardScaler and QuantileTransformer outperformed the other methods, with testing average metric values ranging from 67.8% to 73.7% across different models and oversampling approaches. In contrast, NoScaler and RobustScaler had the lowest performance, with metric values ranging from 64.8% to 73.5%.

For the log-scaled variables, RF and GB were the top-performing classifiers, with average metric values ranging from 66.9% to 74.7% across different oversampling and standardization approaches. In contrast, SVM and QDA had the worst performance, with metrics ranging from 62.4% to 74.2%. Among the oversampling methods, RandomOverSampler and SMOTE performed best, with metrics between 66.7% and 73.1%, while NoOversampling and KMeansSMOTE showed the lowest performance, with values between 65.5% and 73.7%. Regarding standardization techniques, RobustScaler and StandardScaler, as with the non-log-scaled variables, were the best, achieving metric values between 67.3%. and 73.5%. The worst-performing methods were NoScaler and QuantileTransformer, with metric values ranging from 63.9% to 72.6%.

The best-performing combination for both log-scaled and non-log-scaled variables was the RF model using either the RandomOverSampler oversampling using (1) no standardization, (2) StandardScaler standardization. It achieved F1-scores of 76.0% (SD = 1.3%) and 75.9% (SD = 1.9%) for non-log-scaled and log-scaled variables, respectively. Since the standardization method that presented better overall results was the StandardScaler, we will use the model combination of RF with StandardScaler and RandomOverSampler method for further analysis as the best model combination. The worst-performing combination was the SVM model not using both oversampling and standardization, with F1-scores of 37.6% (SD = 2.7%) and 38.8% (SD = 3.5%) for non-log-scaled and log-scaled variables, respectively. Log-scaled variables presented overall better results, having an averaged difference of 0.2% on the F1-scores when compared to the non-log-scaled variables.

The non-log-scaled model showed that the most common misclassifications were: (1) Velociraptor lateral teeth as Saurornitholestes lateral teeth, Suchomimus lateral teeth as Baryonyx lateral teeth, (2) Tyrannosaurus mesial teeth as Tyrannosaurus lateral teeth, (3) Sinraptor teeth as Ceratosaurus lateral teeth or as Allosaurus mesial, (4) Monolophosaurus as Alioramus lateral or as Allosaurus mesial, (5) Masiakasaurus lateral as Saurornitholestes lateral, (6) Majungasaurus mesial as Allosaurus mesial or as Majungasaurus lateral , (7) Gorgosaurus lateral as Allosaurus lateral, (8) Fukuiraptor lateral as Alioramus lateral or as Monolophosaurus, (9) Dromaeosaurus lateral as Saurornitholestes lateral, (10) Daspletosaurus lateral as Acrocanthosaurus lateral or as Tyrannosaurus lateral, (11) Coelophysis lateral as Richardoestesia lateral, (12) Ceratosaurus lateral as Allosaurus lateral, (13) Australovenator as Allosaurus lateral, (14) Allosaurus lateral as Allosaurus mesial, (15) Allosaurus mesial as Allosaurus lateral, (16) Alioramus lateral as Fukuiraptor lateral or as Allosaurus lateral, (17) Albertosaurus lateral as Gorgosaurus lateral (Fig. 6).

Even-though the confusion matrix only shows the predicted class, the three most probable classes were evaluated in the cases were the prediction was incorrect for the whole dataset, resulting in a total of 98.5% correct predictions.

The most influential variable for the predictions was CBW, with a mean importance of 0.22 (Fig. 7), followed by DC, with a mean importance of 0.20 (Fig. 7).

The first branches of the decision trees show that the first criteria to do the classification was if DC is lower or equal to 0.675, the second decision was if (1) CBL was lower or equal to −0.608, (2) DC lower or equal to −0.424 (Fig. 8).

The best performing model was used to classify the set of isolated theropod teeth from Guimarota fossil site and all of the teeth were attributed to their corresponding initial classification in the three most probable classifications, except Ceratosaurus lateral (Table 5).

Classification at higher taxonomic level

The dataset initially included 1,371 theropod teeth from 22 taxa, but after processing, it contained 1,119 isolated theropod teeth belonging to twenty-three combinations of taxa and tooth positions. Two groups were predominant: Tyrannosauridae lateral (175 observations) and Dromaeosauridae lateral (295 observations) (Fig. 9). While the best-performing combination for the non-log-scaled variables was the RF model using StandardScaler and RandomOverSampler technique, achieving F1-scores of 72.0%, the best performing model for log-scaled variables was the RF using QuantileTransformer and RandomOverSampler technique, achieving F1-scores of 72.1%. In contrast, the worst-performing combination for the non-log-scaled variables was the SVM model with KMeansSMOTE oversampling and no standardization, with F1-scores of 44.7%, and for the non-log-scaled variables it was the QDA using BorderlineSMOTE oversampling and no standardization, achieving values of 47.8%. As previously seen on the genus-level classification, the model performance varied significantly across oversampling and scaling techniques. RF and GB were again the best performing models, with consistent results between different oversampling and standardization techniques. GB benefited from RandomOversampler and both RF and GB benefited from QuantileTransformer. NN and SVM improved their performance by using standardization methods (e.g., StandardScaler and RobustScaler), but they remained less effective when compared to RF and GB. QDA and SVM presented the worst metric results and an increase in the metric results when using standardization techniques, namely with StandardScaler. (see Fig. 10 for F1-score and the ShinyApp for the other metrics results). The overall difference for the F1-scores between the log-scaled and non-log-scaled datasets was −0.5%. We will use for further analysis the log-scaled RF model using the QuantileTransformer standardization and RandomOverSampler, since it was the combination that presented best metric results.

The log-scaled model showed that the most common misclassifications were: (1) Noasauridae lateral teeth as Dromaeosauridae lateral, (2) Tyrannosauroidea mesial as Tyrannosauroidea lateral, (3) Carcharodontosauridae mesial as Carcharodontosauridae lateral or as Tyrannosauroidea lateral, (4) Non-abelisauroid Ceratosauria lateral as Metriacanthosauridae lateral or as Tyrannosauroidea lateral, (5) Abelisauridae mesial as Abelisauridae lateral or as Allosauridae lateral, (6) Metriacanthosauridae lateral as non-abelisauroid Ceratosauria lateral, (7) Dromaeosauridae mesial as Dromaeosauridae lateral, (10) Non-averostran Neotheropoda lateral as Dromaeosauridae lateral, (11) Non-abelisauroid Ceratosauria mesial as Allosauridae lateral, (12) Neovenatoridae lateral as Non-spinosaurid Megalosauroidea lateral or as Dromaeosauridae lateral, (13) Basalmost Theropoda as Dromaeosauridae lateral, (14) Abelisauridae mesial as Abelisauridae lateral or as Allosauridae mesial (Fig. 11).

Even-though the confusion matrix only shows the predicted class, the three most probable classes were evaluated in the cases were the prediction was incorrect for the whole dataset, resulting in a total of 97.9% correct predictions.

Similar to the genus-level models, the most important variable for these predictions was CBR, with a mean importance of 0.14 (Fig. 12), followed by DC with a mean importance of 0.12 (Fig. 12).

The first branches of the decision trees show that the first criteria to do the classification was if log-scaled CH is lower or equal to −0.218, the second decision was if (1) log-scaled CBR was lower or equal to 1.477, (2) log-scaled DC lower or equal to 0.789 (Fig. 13).

The best performing model was used to classify the set of isolated theropod teeth from Guimarota fossil site and all of the teeth were attributed to their corresponding classification in the three most probable classifications, except Tyrannosauridae mesial (Table 6)

Classification algorithms

In this study, RF and GB were the most robust models, presenting the best averaged performance and consistency across genus and higher taxonomic-level classifications regardless of oversampling- log-scaling or standardization. Their tree-based structure, that partitions data through feature thresholds rather than distance metrics, allows a more robust and consistent performance of these models, even when data scale changes. StandardScaler presented better averaged metrics for most of the classification tasks, either when using or not log-scaled variables and for the classification at both genus and higher taxonomic levels. The results obtained show that this method can effectively handle outliers and non-normal data distributions, since it standardizes features by removing the mean and scaling to unit variance.

Models like SVM and NN presented more variability, as their performance depends on standardized inputs to avoid skewed distance metrics, benefiting significantly from combined standardization and oversampling. An example of the influence of standardization and oversampling on the model results is the SVM for the genus classifier using the non-log-scaled variables, where the F1- increases 36.7%, from 37.6% (not using any oversampling or standardization technique) to 74.3% (using the BorderlineSMOTE method and the RobustScaler).

Despite the goal of addressing class imbalance, oversampling provided limited improvements, suggesting potential overfitting, where synthetic samples may overly tailor models to the training set without generalizing well to the test data. Tree-based models easily adapted to different scenarios, presenting good performance results without the need of oversampling due to their natural handling of class imbalance through ensemble structures. One of the possible reasons why oversampling did not increase most of the models’ performance is that the oversampled class will be less variable than a class that consists entirely of real data, increasing the possibility of the model not being able to generalize to new data. Standardization played a key role in improving distance-sensitive models like SVM, which relies on balanced feature scales to establish unbiased class boundaries effectively. These findings highlight that while oversampling can benefit SVM and similar models, tree-based models, such as RF and GB, perform optimally across variable scales without oversampling, suggesting that these latter models would be the preferred classifiers in imbalanced datasets.

While using the permutation importance, the two variables with more weight for the predictions in the case of the datasets used in this work were: DC and CBW for the genus level classification; CBR and DC for the higher taxonomic level classification. These findings suggest that variations in denticle density can help differentiate between several genus and clades, making it a valuable trait for classification of isolated theropod teeth as was previously highlighted (Hendrickx et al., 2019). Since CBR is a measure of crown shape proportions, it could capture broader evolutionary trends and ecological adaptations, which are shared among related higher taxonomic level groups. Therefore, theropod higher taxonomic levels with specific dietary niches, such as the robust teeth of apex predators versus the narrow, conical teeth of piscivorous taxa, might show distinct patterns in CBR that differentiate groups like Tyrannosauroidea and Spinosauridae. On the other hand, CBW might have provided finer distinctions within higher taxonomic levels, considering subtle variations in tooth base robustness related to feeding mechanics or dietary specialization unique to specific genera.

In the results section it was possible to see quite common misclassifications between some theropod clades such as Noasauridae lateral and Dromaeosauridae lateral or Allosauridae mesial and Allosauridae lateral (Fig. 14). This misclassification often happens due to similarities in the overall teeth morphology of these theropod clades. Hendrickx, Mateus & Araújo (2014) discusses that Megalosauridae, Ceratosauridae, Abelisauridae, Allosauridae, and Neovenatoridae show considerable overlap, and that significant overlapping also exists for Tyrannosauridae and Carcharodontosauridae, making the distinction between teeth of megalosaurids from the teeth of other similarly sized theropods particularly difficult. The authors also discussed the presence of overlap among different allosauroid taxa, including Allosaurus and Acrocanthosaurus, showing that the distinction between mesial and lateral teeth is not that clear in this clade. Hendrickx, Mateus & Araújo (2014) also divided theropod teeth into four morphospace areas: (1) taxa with small teeth and large denticles (Troodontidae); (2) taxa bearing relatively small teeth and small denticles (non-neotheropod Theropoda, Coelophysoidea, Noasauridae, and Dromaeosauridae); (3) taxa possessing large teeth and minute denticles (Spinosauridae); and (4) ziphodont taxa having relatively large teeth and large denticles (non-noasaurid Ceratosauria, Megalosauridae, Allosauroidea, and Tyrannosauridae). The study argues that overlapping exists between each of these areas, and clades bearing small teeth/denticles and large teeth/ denticles show considerable overlap.

It is important to note that only a small subset of the available variables were used in this analysis, since we only used the variables with less than 15% missing values. While this ensured that the data used was consistent and widely applicable across most teeth, it introduced a possible limitation. Variables that can be critical for identifying certain genus but are incomplete or difficult to measure in other genera are excluded. This restriction could lead to an increased difficulty of the methods to correctly identify genera that are dependent on other variables. This increases the need for researchers to evaluate the reasonability of the predictions, and the importance of evaluating the three most probable classes instead of looking just for the most probable prediction.

Application of ML classifiers to new findings of isolated theropod teeth

The results obtained from our study can be directly applied to the analysis of new findings of isolated theropod teeth, using the trained models and standardizations, or using the same methodologies but training the models from scratch. On the other hand, it is important to note that the training contains a limited number of observations per genus, making it difficult for the models to not overfit to the dataset. Therefore, care should be taken if these models are applied to unseen data from different populations. The following guidelines aim to assist researchers in effectively using the techniques tested in the study of other fossil records.

When classifying isolated theropod teeth, an initial step involves selecting suitable classification algorithms. Tree-based models, such as RF and GB, are robust options, especially when working with raw datasets where standardization and oversampling are not applied. These algorithms rely on decision trees that split data based on feature thresholds, making them well-suited for non-standardized data since they are not sensitive to feature scaling. When applying pre-trained models, it is essential to replicate the pre-processing steps used during training to maintain consistency. If standardization was applied during training, all features must be standardized in the same way, using identical parameters, to ensure the model’s effectiveness.

When classifying paleontological data, oversampling techniques should be considered if the dataset exhibits class imbalance, a common issue due to over-representation of certain taxa. In our study, oversampling notably improved the performance of SVM models, with SMOTE and RandomOversampler proving most reliable. Specifically, the best model for both genus and higher taxonomic level classification used RandomOversampler. This technique is valuable in balancing datasets where teeth from more abundant taxa, such as Tyrannosaurus or Saurornitholestes in the case of the dataset used in this study, overshadow those of less common genera. Implementing RandomOversampler can thus enhance classification performance analyzes. However, oversampling should be avoided when using pre-trained models exclusively for inference on new data, as altering the class distribution could degrade model accuracy.

The decision to apply log-scaling should be guided by the dataset’s characteristics and the model’s sensitivity to skewed distributions. RF and GB models typically perform similarly on log-scaled and non-log-scaled data, as they are less sensitive to distribution shapes. However, log-scaling can be advantageous for algorithms such as SVM and NN, which may benefit from more symmetrically distributed features. Log-scaling, especially when combined with standardization, can stabilize variance and improve model performance when features vary across several orders of magnitude or exhibit skewed distributions.

Application of the methodology in the sample of isolated theropod teeth from Guimarota fossil site

The four isolated theropod teeth mentioned in the methodology section were used to test the classifiers predictions for unseen isolated theropod teeth. The results of the classification using the best models showed that the first three most probable predictions were similar to the original classification. MG 27808_D118 was attributed to Richardoestesia by Zinke (1998) and the results obtained on our classifier models recovered it as Saurornitholestes and Richardoestesia as the firsts and second highest probabilities (Table 5). The results for the other specimens are somewhat more uncertain, with MG 27805_D193 being classified as Troodon and Saurornitholestes, MG 27782 as Troodon and Dilong, and MG 27768 as Masiakasaurus and Saurornitholestes as the first and second highest probabilities in the model based on the genus-level dataset. Based on the higher-taxonomic level dataset, MG 27805_D193 is classified as Dromaeosauridae and Troodontidae, MG 27782 as Troodontidae and Dromaeosauridae, and MG 27768 as Dromaeosauridae and non-averostran Neotheropoda (Table 6). MG 27805_D193 is associated to a labelled with the identification Ceratosaurus-like. An isolated tooth from Guimarota was described by Rauhut, Martin & Krebs (2000) as belonging to a Ceratosaurus-like theropod. However, this specimen has the distal margin almost straight with the apex positioned at the level or more mesial to the most distal end of the distal carina whereas the crown of MG 27805_D193 is strongly recurved, with the mesial margin strongly convex and the distal one concave. Besides, it has an eight-shaped basal cross-section due to the presence of small depressions centrally positioned in the labial and lingual surfaces, which is a character typical for the lateral dentition of dromaeosaurids. Also the presence of much larger distal than mesial denticles (DSDI > 1.2) is also a feature common in most dromaeosaurids. Therefore, the identification obtained based on the classifier models is compatible with the combination of features shown by the specimen suggesting that it probably belongs to a dromaeosaurid theropod.

The specimens MG 27782 and MG 27768 are associated with labels identifying them to tyrannosaurid theropods. Some isolated teeth attributed to tyrannosauroids were described in the Guimarota fossil site by Zinke (1998) and Rauhut, Martin & Krebs (2000). However, these previously described specimens correspond to teeth from a more mesial position in the tooth row compared to the two here analyzed. Both MG 27782 and MG 27768 show the distal carina strongly displaced to the lingual surface, which is a feature shared with several tyrannosauroid taxa but also with several other theropods, including dromaeosaurids (e.g., Hendrickx, Mateus & Araújo (2015)). Therefore, the taxonomic interpretation of these specimens needs to be reviewed based on more uptodate methodologies. An attribution of these isolated teeth to both non-tyrannosaurid tyrannosauroid and dromaeosaurid theropods as suggested by the results of the analyses here performed cannot be excluded at this moment.