Introduction

The package zoolog includes functions and reference data to generate and manipulate log-ratios (also known as log size index (LSI) values) from measurements obtained on zooarchaeological material. Log ratios are used to compare the relative (rather than the absolute) dimensions of animals from archaeological contexts (Meadow 1999). Essentially, the method compares archaeological measurements to a standard, producing a value that indicates how much larger or smaller the archaeological specimen is compared to that standard. zoolog is also able to seamlessly integrate data and references with heterogeneous nomenclature, which is internally managed by a zoolog thesaurus.

The methods included in the package were first developed in the framework of the ERC-Starting Grant 716298 ZooMWest (PI S. Valenzuela-Lamas), and were first used in the paper (Trentacoste, Nieto-Espinet, and Valenzuela-Lamas 2018). They are based on the techniques proposed by G. G. Simpson (1941) and G. Simpson, Roe, and Lewontin (1960), which calculates log size index (LSI) values as: \[ \mbox{LSI} = \log_{10} x - \log_{10} x_{\text{ref}} = \log_{10}(x/x_\text{ref}), \] where \(x\) is the considered measure value and \(x_\text{ref}\) is the corresponding reference value. Observe that LSI is defined using logarithms with base 10.

Several different sets of standard reference values are included in the package. These standards include several published and widely used biometric datasets (e.g. S. J. Davis (1996), Albarella and Payne (2005)) as well as other less known standards. These references, as well as the data example provided with the package, are based on the measures and measure abbreviations defined in Von den Driesch (1976) and S. Davis (1992).

In general, zooarchaeological datasets are composed of skeletal remains representing many different anatomical body parts. In investigation of animal size, the analysis of measurements from a given anatomical element provides the best control for the variables affecting size and shape and, as such, it is the preferable option. Unfortunately, this approach is not always viable due to low sample sizes in some archaeological assemblages. This problem can be mitigated by calculating the LSI values for measurements with respect to a reference, which provides a means of aggregating biometric information from different body parts. The resulting log ratios can be compared and statistically analysed under reasonable conditions (Albarella 2002). However, length, width, and depth measurements of the anatomical elements still should not be directly aggregated for statistical analysis.

The package includes a zoolog thesaurus to facilitate its usage by research teams across the globe, and working in different languages and with different recording traditions. The thesaurus enables the zoolog package to recognises many different names for taxa and skeletal elements (e.g. “Bos taurus”, “cattle”, “BT”, “bovino”, “bota”). Consequently, there is no need to use a particular, standardised recording code for the names of different taxa or elements.

The package also includes a zoolog taxonomy to facilitate the management of cases recorded as identified at different taxonomic ranks (Species, Genus, Tribe) and their match with the corresponding references.

Installation

You can install the released version of zoolog from CRAN with:

install.packages("zoolog")

And the development version from GitHub with:

install.packages("devtools")
devtools::install_github("josempozo/zoolog@HEAD", build_vignettes = TRUE)

Reference standards

The package zoolog includes several osteometrical references. Currently, the references include reference values for the main domesticates and their agriotypes (Bos, Ovis, Capra, Sus), red deer (Cervus elaphus), Persian fallow deer (Dama mesopotamica), mountain gazelle (Gazella gazella), donkey (Equus asinus), horse (Equus caballus), European rabbit (Oryctolagus cuniculus), and grey wolf (Canis lupus). These are drawn from a variety of publications and resources (see below). In addition, the user can consider other references, or the provided references can be extended and updated integrating newer research data. Submission of extended/improved references is encouraged. Please, contact the maintainer through the provided email address to make the new reference fully accessible within the package.

These references, as well as the data example provided with the package, are based on the measurements and measure abbreviations defined in Von den Driesch (1976) and S. Davis (1992). Please, note that Davis’ standard SD, equivalent to Von den Driesch’s DD, has been denoted as Davis.SD in order to resolve its incompatibility with Von den Driesch’s SD. This affects only Davis’ sheep reference and was introduced in Release 1.0.0.

The predefined reference sets included in zoolog are provided in the named list reference, currently comprising the following 4 sets:

library(zoolog)
str(reference, max.level = 1)
#> List of 4
#>  $ NietoDavisAlbarella:'data.frame': 133 obs. of  4 variables:
#>  $ Basel              :'data.frame': 176 obs. of  4 variables:
#>  $ Combi              :'data.frame': 635 obs. of  4 variables:
#>  $ Groningen          :'data.frame': 246 obs. of  4 variables:

The reference set reference$Combi is the default reference for computing the log ratios, since it includes the most comprehensive reference for each species.

The package also includes a referencesDatabase collecting the taxon-specific reference standards from all the considered resources. Each reference set is composed of a different combination of taxon-specific standards selected from this referencesDatabase. The selection is defined by the data frame referenceSets:

A detailed description of the reference data, including its structure, properties, and considered resources can be found in the ReferencesDatabase help page.

Thesaurus

A thesaurus set is defined in order to make the package compatible with the different recording conventions and languages used by authors of zooarchaeological datasets. This enables the function LogRatios to match values in the user’s dataset with the corresponding ones in the reference standard, regardless of differences in nomenclature or naming conventions, as long as both terms are included in the relevant thesaurus. The thesaurus also allows the user to standardize the nomenclature of the dataset if desired.

The user can also use other thesaurus sets or modify the provided one. In this latter case, we encourage the user to contact the maintainer at the provided email address so that the additions can be incorporated into the new versions of the package.

Currently, the zoolog thesaurus set includes four thesauri:

identifierThesaurus

For the column names that identify the variables used in computing the log ratios. It includes the categories Taxon, Element, Measure, and Standard. Each category provides a series of equivalent names. For instance, Taxon includes the options Taxon, TAX, species, animal, Specie, GenusSpecies, Especie, Espece. This thesaurus is case, accent, and punctuation insensitive, so that, for instance, “Especie” is equivalent to “ESPECIE” or “Espècie”.

taxonThesaurus

For the names of the different taxa when recording animal bones. The current categories are Bos taurus, Bos primigenius, Bos, Ovis aries, Ovis orientalis, Ovis, Capra hircus, Capra aegagrus, Capra, Caprini, Sus domesticus, Sus scrofa, Sus, Cervus elaphus, Cervus, Dama mesopotamica, Dama, Gazella gazella, Gazella, Equus asinus, Equus caballus, Equus, Oryctolagus cuniculus, Oryctolagus, Canis familiaris, Canis lupus, Canis, each with different equivalent names. This thesaurus is case, accent, and punctuation insensitive, so that, for instance, “Bos” is equivalent to “bos”, “Bos.”, or “Bos ”.

elementThesaurus

Names of anatomical elements when recording animal bones. It currently includes 53 categories (scapula, humerus, radius, …), each with different equivalent names. This thesaurus is case, accent, and punctuation insensitive.

measureThesaurus

Names of the measurements. While the English abbreviations from Von den Driesch (1976) and S. Davis (1992) are widely used in published literature, this thesaurus enables other nomenclatures (e.g. original German abbreviations in Von den Driesch (1976)) to be included. This thesaurus is case sensitive.

Taxonomy

A taxonomy for the most typical zooarchaeological taxa has been introduced in the zoolog major release 1.0.0. The taxonomic hierarchy is structured up to the family level:

This taxonomy is intended to facilitate the management of cases recorded as identified at different taxonomic ranks and their match with the corresponding references.

The taxonomy is integrated in the function LogRatios, enabling it to automatically match different species in data and reference that are under the same genus. For instance, data of Bos taurus can be matched with reference of Bos primigenius, since both are Bos. It also enables the function LogRatios to detect when a case taxon has been only partially identified and recorded at a higher taxonomic rank, such as tribe or family, and to suggest the user the set of possible reference species.

Besides, it is complemented with a series of functions enabling the user to query for the subtaxonomy or the set of species under a queried taxon at any taxonomic rank.

Functions

The full list of functions is available under the zoolog help page. We list them here sorted by their prominence for a typical user, and grouped by functionality:

LogRatios

It computes the log ratios of the measurements in a dataset relative to standard reference values. By default reference$Combi is used. The function includes the option ‘joinCategories’ allowing several taxa (typically Ovis, Capra, and unknown Ovis/Capra) to be considered together with the same reference taxon.

Note that without using ‘joinCategories’ any taxa not part of the selected reference set will be excluded. For instance, if using reference$NietoDavisAlbarella, log ratios for goats will not be calculated unless ‘joinCategories’ is set to indicate that the Ovis aries standard should also be applied to goats.

CondenseLogs

It condenses the calculated log ratio values into a reduced number of features by grouping several measure log ratios and selecting or calculating a representative feature value. By default the selected groups represent a single dimension, i.e. Length, Width, and Depth. Only one feature is extracted per group. Currently, two methods are possible: “priority” (default) or “average”.

This operation is motivated by two circumstances. First, not all measurements are available for every bone specimen, which obstructs their direct comparison and statistical analysis. Second, several measurements can be strongly correlated (e.g. SD and Bd both represent bone width). Thus, considering them as independent would produce an over-representation of bone remains with multiple measurements per axis. Condensing each group of measurements into a single feature (e.g. one measure per axis) alleviates both problems.

The default method (“priority”), selects the first available log ratio in each group. Besides, CondenseLogs employs the following by-default group and prioritization introduced in Trentacoste, Nieto-Espinet, and Valenzuela-Lamas (2018): Length considers in order of priority GL, GLl, GLm, and HTC. Width considers in order of priority BT, Bd, Bp, SD, Bfd, and Bfp. Depth considers in order of priority Dd, DD, BG, and Dp. This order maximises the robustness and reliability of the measurements, as priority is given to the most abundant, more replicable, and less age dependent measurements. But users can set their own features with any group of measures and priorities. The method “average” extracts the mean per group, ignoring the non-available log ratios.

RemoveNACases

It removes the cases (table rows) for which all measurements of interest are non-available (NA). A particular list of measurement names can be explicitly provided or selected by a common initial pattern (e.g. prefix). The default setting removes the rows with no available log ratios to facilitate subsequent analysis of the data.

InCategory

It checks if an element belongs to a category according to a thesaurus. It is similar to base::is.element, returning a logical vector indicating if each element in a given vector is included in a given set. But InCategory checks for equality assuming the equivalencies defined in the given thesaurus. It is intended for the user to easily select a subset of data without having to standardize the analysed dataset.

Nomenclature standardization

This includes two functions enabling the user to map data with heterogeneous nomenclature into a standard one as defined in a thesaurus:

• StandardizeNomenclature standardizes a character vector according to a given thesaurus.
• StandardizeDataSet standardizes column names and values of a data frame according to a thesaurus set.

Subtaxonomy

This includes three functions enabling the user to query for some information on the subtaxonomy under a queried taxon at any taxonomic rank:

• Subtaxonomy provides the subtaxonomy dataframe collecting the species (rows) included in the queried taxon, and the taxonomic ranks (columns) up to its level.
• SubtaxonomySet provides the set (character vector without repetions) of taxa, at any taxonomic rank, under the queried taxon.
• GetSpeciesIn provides the set of species included in the queried taxon.

By default, the subtaxonomy information is extracted from the zoolog taxonomy.

AssembleReference

It allows the user to build new references assembling the desired taxon-specific references included in the zoolog referencesDataSet or in any other provided by the user.

Thesaurus readers and writers

This includes functions to read and write a single thesaurus (ReadThesaurus and WriteThesaurus) and a thesaurus set (ReadThesaurusSet and WriteThesaurusSet).

Thesaurus management

This includes functions to modify and check thesauri:

• NewThesaurus generates an empty thesarus.
• AddToThesaurus adds new names and categories to an existing thesaurus.
• RemoveRepeatedNames cleans a thesarus from any repeated names on any category.
• ThesaurusAmbiguity checks if there are names included in more than one category in a thesaurus.

Examples

The following examples are designed to be read and run sequentially. They represent a possible pipeline, meaningful for the processing and analysis of a dataset. Only occasionally, a small diversion is included to illustrate some alternatives.

library(zoolog)
dataFile <- system.file("extdata", "dataValenzuelaLamas2008.csv.gz", 
                        package = "zoolog")
data = read.csv2(dataFile,
                 quote = "\"", header = TRUE, na.strings = "",
                 encoding = "UTF-8")
knitr::kable(head(data)[, -c(6:20,32:64)])

dataWithLog <- LogRatios(data)
#> Warning: Reference for Sus scrofa used for cases of Sus domesticus.
#>    Reference for Sus scrofa used for cases of Sus.
#>    Reference for Oryctolagus cuniculus used for cases of Oryctolagus.
#>    Reference for Canis lupus used for cases of Canis.
#>    Reference for Ovis aries used for cases of Ovis.
#>    Set useGenusIfUnambiguous to FALSE if this behaviour is not desired.
#> Warning: Data includes some cases recorded as
#>     * Equus (which is a Genus)
#>       for which the reference for Equus asinus or Equus caballus could be used.
#>     * Caprini (which is a Tribe)
#>       for which the reference for Ovis aries or Capra hircus could be used.
#>    Set joinCategories as appropriate if you want to use any of them.
knitr::kable(head(dataWithLog)[, -c(6:20,32:64)])

AScases <- InCategory(dataWithLog$Os, "astragalus", zoologThesaurus$element)
knitr::kable(head(dataWithLog[AScases, -c(6:20,32:64)]))

GLVariants <- list(c("GL", "GLl", "GLpe"))
dataWithLog <- LogRatios(data, mergedMeasures = GLVariants)
knitr::kable(head(dataWithLog[AScases, -c(6:20,32:64)]))

caprineCategory <- list(ovar = SubtaxonomySet("caprine"))
dataWithLog <- LogRatios(data, joinCategories = caprineCategory, mergedMeasures = GLVariants)
knitr::kable(head(dataWithLog)[, -c(6:20,32:64)])

SubtaxonomySet("caprine")
#> [1] "Ovis aries"      "Ovis orientalis" "Capra hircus"    "Capra aegagrus" 
#> [5] "Ovis"            "Capra"           "Caprini"

dataWithLogPruned=RemoveNACases(dataWithLog)
knitr::kable(head(dataWithLogPruned[, -c(6:20,32:64)]))

write.csv2(dataWithLogPruned, "myDataWithLogValues.csv", 
           quote=FALSE, row.names=FALSE, na="", 
           fileEncoding="UTF-8")

dataWithSummary <- CondenseLogs(dataWithLogPruned)
knitr::kable(head(dataWithSummary)[, -c(6:20,32:64,72:86)])

dataStandardized <- StandardizeDataSet(dataWithSummary)
knitr::kable(head(dataStandardized)[, -c(6:20,32:64,72:86)])

dataOC <- subset(dataWithSummary, InCategory(Especie, 
                                             SubtaxonomySet("caprine"),
                                             zoologThesaurus$taxon))
knitr::kable(head(dataOC)[, -c(6:20,32:64)])

dataOCStandardized <- StandardizeDataSet(dataOC)
knitr::kable(head(dataOCStandardized)[, -c(6:20,32:64)])

dataOCWithWidth <- RemoveNACases(dataOCStandardized, measureNames = "Width")
dataOCWithLength <- RemoveNACases(dataOCStandardized, measureNames = "Length")

ggplot(dataOCStandardized, aes(x = Site, y = Width)) +
  geom_boxplot(outlier.shape = NA, na.rm = TRUE) +
  geom_jitter(width = 0.2, height = 0, alpha = 1/2, color = 4, na.rm = TRUE) +
  theme_bw() +
  ggtitle("Caprine widths") +
  ylab("Width log-ratio") +
  coord_flip()

ggplot(dataOCStandardized, aes(x = Site, y = Length)) +
  geom_boxplot(outlier.shape = NA, na.rm = TRUE) +
  geom_jitter(width = 0.2, height = 0, alpha = 1/2, color = 4, na.rm = TRUE) +
  theme_bw() +
  ggtitle("Caprine lengths") +
  ylab("Length log-ratio") +
  coord_flip()

ggplot(dataOCStandardized, aes(Width)) +
  geom_histogram(bins = 30, na.rm = TRUE) +
  ggtitle("Caprine widths") +
  xlab("Width log-ratio") +
  facet_grid(Site ~.) +
  theme_bw() +
  theme(panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5, size = 14),
        axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10),
        axis.text = element_text(size = 10) ) +
  scale_y_continuous(breaks = c(0, 10, 20, 30))

levels0 <- unique(dataOCStandardized$Taxon)
levels0
#> [1] "Ovis aries"   "Capra hircus" "Caprini"

dataOCStandardized$Taxon <- factor(dataOCStandardized$Taxon, 
                                   levels = levels0[c(1,3,2)])
levels(dataOCStandardized$Taxon)
#> [1] "Ovis aries"   "Caprini"      "Capra hircus"

Ocolour <- c("#A2A475", "#D8B70A", "#81A88D")
ggplot(dataOCStandardized, aes(x=Site, y=Width)) + 
  geom_boxplot(aes(fill=Taxon), 
               notch = TRUE, alpha = 0, lwd = 0.377, outlier.alpha = 0,
               width = 0.5, na.rm = TRUE,
               position = position_dodge(0.75),
               show.legend = FALSE) + 
  geom_point(aes(colour = Taxon, shape = Taxon), 
             alpha = 0.7, size = 0.8, 
             position = position_jitterdodge(jitter.width = 0.3),
             na.rm = TRUE) +
  scale_colour_manual(values=Ocolour) +
  scale_shape_manual(values=c(15, 18, 16)) +
  theme_bw(base_size = 8) +
  ylab("LSI value") +
  ggtitle("Sheep/goat LSI width values") 

TaxonSiteWidthHist <- ggplot(dataOCStandardized, aes(Width, fill = Taxon)) + 
  geom_histogram(bins = 30, alpha = 0.5, position = "identity") + 
  ggtitle("Sheep/goat Widths") + facet_grid(Site ~ Taxon)
TaxonSiteWidthHist

TaxonSiteWidthHist <- ggplot(dataOCStandardized, aes(Width, fill = Taxon)) + 
  geom_histogram(bins = 30, alpha = 0.5, position = "identity") + 
  ggtitle("Sheep/goat Widths") + facet_grid(~Site)
TaxonSiteWidthHist

t.test(Length ~ Site, dataOCStandardized, 
       subset = Site %in% c("OLD", "ALP"),
       na.action = "na.omit")
#> 
#>  Welch Two Sample t-test
#> 
#> data:  Length by Site
#> t = -12.23, df = 4.2233, p-value = 0.0001865
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -0.11262797 -0.07165242
#> sample estimates:
#> mean in group ALP mean in group OLD 
#>       -0.04589475        0.04624545

t.test(Width ~ Site, dataOCStandardized, 
       subset = Site %in% c("OLD", "ALP"),
       na.action = "na.omit")
#> 
#>  Welch Two Sample t-test
#> 
#> data:  Width by Site
#> t = -3.0613, df = 8.8203, p-value = 0.01386
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -0.09749890 -0.01448767
#> sample estimates:
#> mean in group ALP mean in group OLD 
#>       -0.03922705        0.01676624

library(stats)
pairwise.t.test(dataOCStandardized$Width, dataOCStandardized$Site, 
                pool.sd = FALSE)
#> 
#>  Pairwise comparisons using t tests with non-pooled SD 
#> 
#> data:  dataOCStandardized$Width and dataOCStandardized$Site 
#> 
#>     ALP   OLD  
#> OLD 0.042 -    
#> TFC 0.450 0.051
#> 
#> P value adjustment method: holm

How to cite the package zoolog

We have invested a lot of time and effort in creating zoolog. Please cite the package if you publish an analysis or results obtained using zoolog. For example, “Log ratios calculation and analysis was done using the package zoolog (Pozo et al. 2022) in R 4.0.3 (R Core Team 2020).”

To get the details of the most recent version of the package, you can use the R citation function: citation("zoolog").

When publishing it is also recommended that the references for standards are properly cited, as well as any details on selecting feature values. The details on the source of each of the reference standards are given in the help page referencesDatabase. For instance, if using the zoolog reference$NietoDavisAlbarella and the default selection method of features, a fair description may be: “Published references for cattle (Nieto-Espinet 2018), sheep/goat (S. J. Davis 1996) and pigs (Albarella and Payne 2005) were used as standards. One length and one width log ratio value from each specimen were included in the analysis, with values selected following the default zoolog ‘priority’ method (Trentacoste, Nieto-Espinet, and Valenzuela-Lamas 2018): length values - GL, GLl, GLm, HTC; width values - BT, Bd, Bp, SD, Bfd, Bfp.”