20 Visualisation of lead isotope data

Authors

Yiu-Kang Hsu

Published

October 5, 2026

20.1 Learning objective

By the end of this unit you are aware about the benefits and limitations of the different ways lead isotope data can be visualised.

20.2 Prior knowledge

This section assumes familiarity with the prior learning materials on Pb isotope geochemistry, esp. Chapter 3.

20.3 Material

This section discusses different types of plots. Interactive examples of these plots allows you to explore their suitability for the different research questions on your own.

20.4 Learning content

Because lead isotopes are represented by three independent ratios, e.g., ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb, they can be visualised in a three dimensional geometric space. However, often only two dimensions can be represented in one plot. In addition, different ways of optically grouping and highlighting groups within the Pb isotope space exist. It is important to keep in mind that every plot has its pros and cons and there is no general consensus which is the best presentation.

All plots were created in R with support of different packages. The source of the respective chunks is provided in collapsible sections. The entire source can be copied from the raw version of this chapter through the code tools next to the chapter heading.

Load R packages and datasets

library(ggplot2)     # creating static plots
library(plotly)      # creating interactive plots
library(ggtern)      # ternary plots
library(ks)          # calculation of kernel densities
library(rgl)         # creating interactive 3D plots
library(viridisLite) # colour blind-friendly colour schemes

setupKnitr(autoprint = TRUE)

The same subset of a real archaeological dataset will be used as example to demonstrate how this data can be treated in different ways. The dataset is available in the GitHub repository of this book and was originally published by (Bode 2008). To avoid overloading of the plots, the subset is restricted to the archaeological objects and reference data from the region “Northern Eifel”, “Sauerland”, and “Bergisches Land”, all located in the West of Germany.

Import and filter data

Roman_case <- read.csv("example_dataset/Roman.csv", encoding="UTF-8", header = TRUE) |>
  subset(Region %in% c("Northern Eifel", "Sauerland", "Bergisches Land", "Roman object"))

In addition, the plots share several stylistic elements, such as the colour scheme. Therefore, they can be defined for all plots.

Define common style elements

# define plotting order of regions
Roman_case$Region <- factor(Roman_case$Region, levels = c('Northern Eifel', 'Sauerland', 'Bergisches Land', 'Roman object')) 

# define colours 
colours <- c(viridis(3), "black")|>
  setNames(c('Northern Eifel', 'Sauerland', 'Bergisches Land', 'Roman object'))

# define symbols
# in R, symbols are defined by numbers, see e.g., https://cran.r-project.org/web/packages/ggplot2/vignettes/ggplot2-specs.html#point
symbols <- c(23,24,22,3) |>
  setNames(c('Northern Eifel', 'Sauerland', 'Bergisches Land', 'Roman object'))

20.4.1 Binary scatter plot

The bi-plot (Figure 20.1) is by far the most common option to display lead isotope data. Since there are four isotopes of Pb, twelve combinations of isotopic ratios can be derived. The use of paired ratios depends on the instruments used and the scientific disciplines of the studies. In the early days, Pb isotopic ratios were often reported based on ²⁰⁶Pb-based ratios as ²⁰⁴Pb could not be measured precisely. However, in the 2000s, the advent of the multi-collector mass spectrometer (MC-ICP-MS) and the double- or triple-spiked technique created a huge amount of Pb isotope data with precisely measured ²⁰⁴Pb. Conventionally, environmental science tends to use the ratios based on ²⁰⁶Pb, which however generates plots with linear patterns and thus a low discrimination power (Ellam 2010). In geological literature, ratios based on ²⁰⁴Pb are commonplace which enable a better visualisation of system closure time (or model age) and U-Th-Pb composition (or µ and κ) of parental source(s) (Albarède et al. 2012). However, it has to be kept in mind that all two-dimensional plots incompletely represent a dataset. All twelve combination plots are suggested to be tested to view the full isotopic extent of ore deposits (Albarède et al. 2020). Ideally, the Pb isotopic ratios should be considered in a three-dimensional space.

Show the code

# select the region and all columns with isotope ratios (= all columns starting with Pb followed by three characters and a period as replacement character for the slash)  
data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}\\.", names(Roman_case), value = TRUE)))

colnames(data) <- c('Region', '204Pb/206Pb', '207Pb/206Pb', '208Pb/206Pb', '206Pb/204Pb', '207Pb/204Pb', '208Pb/204Pb')

# create plot
plt <- plot_ly(data = data, x = ~`206Pb/204Pb`, y = ~`207Pb/204Pb`, 
        split = ~Region,  color = ~Region, symbol = ~Region, text = ~Region,
        colors = colours, alpha = 0.75, symbols = symbols, 
        marker = list(line = list(color = 'rgb(0, 0, 0)', width = 1)),
        hovertemplate = paste('Region: %{text}<br>',
                              'x: %{x}<br>',
                              'y: %{y}', 
                              '<extra></extra>')
        )

# set-up drop down menus for the x and y variables 

  # function to create the drop down menu elements
  # adapted from a function provided on https://stackoverflow.com/questions/75596474/add-a-legend-to-a-plotly-plot-with-drop-down-menus-in-r
  btn <- function(xLoc, m, data, xOry) { # x btn position, method, names, x or y
    list(type = "list", x = xLoc, xanchor = "left", y = 1.2, 
         buttons = lapply(
           names(data[[1]])[-1], function(k) {           # iterate over names
             switch (xOry,
                     x = {args <- list(    # a list for each trace (each legend entry)
                       list(x = list(data[[1]][, k], data[[1]][, k], data[[1]][, k], data[[1]][, k])),  # 4 lists, 4 traces
                       list(xaxis = list(title = k)))
                     }, 
                     y = {args = list(    # a list for each trace (each legend 
                       list(x = list(data[[1]][, k], data[[1]][, k], data[[1]][, k], data[[1]][, k])),  # 4 lists, 4 traces
                       list(yaxis = list(title = list(text = k))))}, 
                     stop("xOrY must be either 'x' or 'y'.")
             )
             list(method = m, label = k, args = args)
           }
         )
    )
  }

  # split regions into separate datasets 
  regions <- split(data, data$Region) 

# add dropdown menus to plot
plt <- plt %>% layout(
  xaxis = list(title = "206Pb/204Pb"),
  yaxis = list(title = "207Pb/204Pb"),
  updatemenus = list(btn(.2, "update", regions, "x"),
                     btn(.8, "update", regions, "y")), 
  showlegend = TRUE,
  annotations = list(
    list(
      text = "<b>X-Axis:</b>", x=0.04, y=1.18, 
      xref='paper', yref='paper',xanchor = "left", showarrow=FALSE
    ),
    list(
      text = "<b>Y-Axis:</b>", x=0.63, y=1.18, 
      xref='paper', yref='paper',xanchor = "left", showarrow=FALSE
    )
  )
)

# draw plot
plt

Figure 20.1: A binary plot of Roman lead artefacts in comparison with ore districts in Germany.

20.4.2 Bi-plot using geological-informed parameters

Instead of isotopic ratios, Albarède et al. (2012) advocate the use of calculated geological model parameters, namely the model age (T), U/Pb (μ), and Th/U (κ) to discriminate potential ore sources in provenance studies (Figure 20.2). As shown in chapter 3, ²⁰⁶Pb, ²⁰⁷Pb, and ²⁰⁸Pb are generated by radioactive decay of their parental isotopes ²³⁸U, ²³⁵U, and ²³²Th, respectively. We can therefore calculate the model age, ²³⁸U/²⁰⁴Pb and ²³²Th/²³⁸U from the Pb isotope ratios determined for a given sample using the equations provided in Albarède et al. (2012) or any other of the Pb isotope models mentioned in chapter 3 by using, e.g., an R script.

Show the code

# select columns Age, Mu, Kappa, Region and rename them 
data <- subset(Roman_case, select=c('Region', 'Age', 'Mu', 'Kappa')) |>
  na.omit()

colnames(data) <- c('Region', 'Age','U/Pb','Th/U')

# create plot
plt <- plot_ly(data = data, x = ~Age, y = ~`U/Pb`, 
        split = ~Region,  color = ~Region, symbol = ~Region, text = ~Region,
        colors = colours, alpha = 0.75, symbols = symbols, 
        marker = list(line = list(color = 'rgb(0, 0, 0)', width = 1)),
        hovertemplate = paste('Region: %{text}<br>',
                              'x: %{x}<br>',
                              'y: %{y}', 
                              '<extra></extra>')
        )

# set-up drop down menus for the x and y variables. 

  # The function 'btn', defined in the code chunk above, is reused.   

  # split regions into separate datasets 
  regions <- split(data, data$Region) 

# add dropdown menus to plot
plt <- plt %>% layout(
  xaxis = list(title = "Age"),
  yaxis = list(title = "U/Pb"),
  updatemenus = list(btn(.2, "update", regions, "x"),
                     btn(.8, "update", regions, "y")), 
  showlegend = TRUE,
  annotations = list(
    list(
      text = "<b>X-Axis:</b>", x=0.04, y=1.18, 
      xref='paper', yref='paper',xanchor = "left", showarrow=FALSE
    ),
    list(
      text = "<b>Y-Axis:</b>", x=0.63, y=1.18, 
      xref='paper', yref='paper',xanchor = "left", showarrow=FALSE
    )
  )
)

# draw plot
plt

Figure 20.2: A binary plot of Roman lead artefacts using geological parameters of model age (T), U/Pb (μ), and Th/U (κ).

20.4.3 Bi-plot with 90% confidence ellipse

The increasing amount of available Pb isotope analyses resulted in the issue of how to effectively present overlapping data. One way to circumvent this problem is the creation of so-called “ore fields” that represent the extent of Pb isotopic variation given that the data follow a normal distribution. Some researchers have resorted to the 90% confidence ellipse (Figure 20.3) as the reference for ore sourcing (Stos-Gale et al. 1997). This technique was severely criticised by Baxter and Gale (1998), who demonstrated the non-normality of Pb isotope data in a number of instances. The 90% confidence ellipse is no longer considered suitable for representing an ore Pb isotopic population.

Show the code

# select the region and all 204Pb-normalised isotope ratios  
data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}\\.Pb204", names(Roman_case), value = TRUE)))

colnames(data) <- c('Region', '206Pb/204Pb', '207Pb/204Pb', '208Pb/204Pb')

# create plot for 207Pb/204Pb vs 206Pb/204Pb
p_207 <- ggplot(data, aes(x = `206Pb/204Pb`, y = `207Pb/204Pb`, shape = Region, fill = Region)) +
  geom_point(stroke = 0.2,size = 2.5, color = "black", alpha = 0.75) +
  scale_shape_manual(values = symbols) + 
  scale_fill_manual(values = colours) +
  stat_ellipse(level = 0.9, type = "t", linewidth = 0.1, show.legend = NA) +
  labs(x = expression(""^"206"*"Pb/"^"204"*"Pb"), 
       y = expression(""^"207"*"Pb/"^"204"*"Pb")) + 
  theme_bw()

p_207 <- ggplotly(p_207, tooltip = c("x", "y", "fill")) %>% 
  layout(xaxis = list(title = "<sup>206</sup>Pb/<sup>204</sup>Pb"),
         yaxis = list(title = "<sup>207</sup>Pb/<sup>204</sup>Pb"),
         legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175)
         )

# create plot for 208Pb/204Pb vs 206Pb/204Pb
p_208 <- ggplot(data, aes(x = `206Pb/204Pb`, y = `208Pb/204Pb`, shape = Region, fill = Region)) +
  geom_point(stroke = 0.2,size = 2.5, color = "black", alpha = 0.75) +
  scale_shape_manual(values = symbols) + 
  scale_fill_manual(values = colours) +
  stat_ellipse(level = 0.9, type = "t", linewidth = 0.1, show.legend = NA) +
  labs(x = expression(""^"206"*"Pb/"^"204"*"Pb"), 
       y = expression(""^"208"*"Pb/"^"204"*"Pb")) + 
  theme_bw()

p_208 <- ggplotly(p_208, tooltip = c("x", "y", "fill")) %>% 
  layout(xaxis = list(title = "<sup>206</sup>Pb/<sup>204</sup>Pb"),
         yaxis = list(title = "<sup>208</sup>Pb/<sup>204</sup>Pb"),
         legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175)
         )

# draw plots
p_207
p_208

(a) Ellipses based on ²⁰⁶Pb/²⁰⁴Pb vs ²⁰⁷Pb/²⁰⁴Pb.

(b) Ellipses based on ²⁰⁶Pb/²⁰⁴Pb vs ²⁰⁸Pb/²⁰⁴Pb.

Figure 20.3: Binary scatter plots with 90% confidence ellipses around data points.

20.4.4 Binary plot with the kernel density estimation

Due to the non-normality of Pb isotope data, both Baxter et al. (1997) and Scaife et al. (1999) advocated the use of more robust kernel density estimates (KDEs) to display the isotopic extent of orefields (Figure 20.4). KDEs are a non-parametric method to transform continuous data into a smoothed probability density function. KDEs offer three main advantages:

They do not assume the normality of data;
They can produce smoother distributions than conventional histograms, whose appearance is significantly affected by the choices of bin width and the start/end points of bins; and
They can represent data in a multidimensional space and enable users to effectively compare different datasets either graphically or mathematically.

Given these advantages, the KDE method has become popular in recent publications of archaeological sciences (Hsu et al. 2018).

Show the code

# select individual regions/assemblages
regions <- c("Northern Eifel", "Sauerland", "Bergisches Land")

# define levels, to be read as 100% - value, i.e. '5%' indicates 95% contour
levels <- c("5%", "25%", "50%") 

# initiate data
data_kde206v207 <- vector(mode = "list", length = length(regions)) 
names(data_kde206v207) <- regions

data_kde206v208 <- vector(mode = "list", length = length(regions))
names(data_kde206v208) <- regions

for (i in regions) {
  
  # subset data
  data_region <- Roman_case[Roman_case$Region == i, ]
  
  # compute KDE
  data_region_kde206v207 <- ks::kde(data_region[c("Pb206.Pb204", "Pb207.Pb204")],  compute.cont = TRUE,gridsize = 1024)
  data_region_kde206v208 <- ks::kde(data_region[c("Pb206.Pb204", "Pb208.Pb204")],  compute.cont = TRUE,gridsize = 1024)
  
  # compute contour lines
  
  data_region_kde206v207_levels <- data_region_kde206v207$cont[levels]
  
  data_region_kde206v207_contour <- with(data_region_kde206v207, contourLines(x = eval.points[[1]], y = eval.points[[2]], z = estimate, levels = data_region_kde206v207_levels))
  names(data_region_kde206v207_contour) <- paste0(i, seq_along(data_region_kde206v207_contour))
  data_region_kde206v207_contour <- Map(c, data_region_kde206v207_contour, group = names(data_region_kde206v207_contour))
  data_region_kde206v207_contour <- do.call("rbind", lapply(data_region_kde206v207_contour, data.frame))
  data_region_kde206v207_contour$level_percent <- names(data_region_kde206v207_levels)[match(data_region_kde206v207_contour$level, data_region_kde206v207_levels)]
  
  data_region_kde206v208_levels <- data_region_kde206v208$cont[levels]
  
  data_region_kde206v208_contour <- with(data_region_kde206v208, contourLines(x = eval.points[[1]], y = eval.points[[2]], z = estimate, levels = data_region_kde206v208_levels))
  names(data_region_kde206v208_contour) <- paste0(i, seq_along(data_region_kde206v208_contour))
  data_region_kde206v208_contour <- Map(c, data_region_kde206v208_contour, group = names(data_region_kde206v208_contour))
  data_region_kde206v208_contour <- do.call("rbind", lapply(data_region_kde206v208_contour, data.frame))
  data_region_kde206v208_contour$level_percent <- names(data_region_kde206v208_levels)[match(data_region_kde206v208_contour$level, data_region_kde206v208_levels)]
  
  data_kde206v207[[i]] <- data_region_kde206v207_contour
  data_kde206v208[[i]] <- data_region_kde206v208_contour
  
}

# crunch list into data frame
data_kde206v207 <- Map(cbind, data_kde206v207, Region = names(data_kde206v207))
data_kde206v207 <- do.call("rbind", data_kde206v207)

data_kde206v208 <- Map(cbind, data_kde206v208, Region = names(data_kde206v208))
data_kde206v208 <- do.call("rbind", data_kde206v208)

# plots
plot_kde206v207 <- ggplot() + 
  geom_point(data = Roman_case[Roman_case$Region == "Roman object", ], aes(x = Pb206.Pb204, y = Pb207.Pb204, colour = Region), shape = 3, size = 2.3) + 
    geom_polygon(data = data_kde206v207, mapping = aes(x = x, y= y, group = group, colour = Region, fill = Region, alpha = level_percent)) +
  scale_colour_manual(name = NULL, limits = names(colours), values = colours) + 
  scale_fill_manual(name = NULL, limits = names(colours), values = colours) + 
  scale_alpha_discrete(range = c(0.1, 0.3), guide = "none") +
  scale_x_continuous(name = expression(""^"206"*"Pb/"^"204"*"Pb"), limits = c(18, 18.6), breaks = seq(18, 18.6, by = 0.1)) +
  scale_y_continuous(name = expression(""^"207"*"Pb/"^"204"*"Pb"), limits = c(15.5, 15.8), breaks = seq(15.5, 15.8, by = 0.05)) +
  guides(colour = guide_legend(position = "inside")) +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), 
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 14),
        legend.position.inside = c(0.76, 0.11)
        ) 

plot_kde206v208 <- ggplot() + 
    geom_point(data = Roman_case[Roman_case$Region == "Roman object", ], aes(x = Pb206.Pb204, y = Pb208.Pb204, colour = Region), shape = 3, size = 2.3) +
  geom_polygon(data = data_kde206v208, mapping = aes(x = x, y = y, group = group, colour = Region, fill = Region, alpha = level_percent)) +
    geom_point(data = Roman_case[Roman_case$Region == "Roman object", ], aes(x = Pb206.Pb204, y = Pb208.Pb204, colour = Region), shape = 3, size = 2.3) + 
  scale_colour_manual(name = NULL, limits = names(colours), values = colours) + 
  scale_fill_manual(name = NULL, limits = names(colours), values = colours) + 
  scale_alpha_discrete(range = c(0.1, 0.3), guide = "none") +
  scale_x_continuous(name = expression(""^"206"*"Pb/"^"204"*"Pb"), limits = c(18, 18.6), breaks = seq(18, 18.6, by = 0.1)) +
  scale_y_continuous(name = expression(""^"208"*"Pb/"^"204"*"Pb"), limits = c(37.9, 38.8), breaks = seq(38, 39, by = 0.1)) +
  guides(colour = guide_legend(position = "inside")) +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), 
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 14),
        legend.position.inside = c(0.76, 0.11)
        ) 

# convert to plotly
plot_kde206v207 <- ggplotly(plot_kde206v207) %>% 
  layout(xaxis = list(title = "<sup>206</sup>Pb/<sup>204</sup>Pb"),
         yaxis = list(title = "<sup>207</sup>Pb/<sup>204</sup>Pb"),
         legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175)
         )

plot_kde206v208 <- ggplotly(plot_kde206v208)  %>% 
  layout(xaxis = list(title = "<sup>206</sup>Pb/<sup>204</sup>Pb"),
         yaxis = list(title = "<sup>207</sup>Pb/<sup>204</sup>Pb"),
         legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175)
         )

# modify legends
for (i in seq_along(plot_kde206v207$x$data)) {
  if (!grepl("25%|Roman", plot_kde206v207$x$data[[i]]$name)) {plot_kde206v207$x$data[[i]]$showlegend <- FALSE}
  for (k in unique(Roman_case$Region)) {
    if (grepl(k, plot_kde206v207$x$data[[i]]$name)) {
      plot_kde206v207$x$data[[i]]$name <- k
      plot_kde206v207$x$data[[i]]$legendgroup <- k
    }
      
  }
}

for (i in seq_along(plot_kde206v208$x$data)) {
  if (!grepl("25%|Roman", plot_kde206v208$x$data[[i]]$name)) {plot_kde206v208$x$data[[i]]$showlegend <- FALSE}
  for (k in unique(Roman_case$Region)) {
    if (grepl(k, plot_kde206v208$x$data[[i]]$name)) {
      plot_kde206v208$x$data[[i]]$name <- k
      plot_kde206v208$x$data[[i]]$legendgroup <- k
    }
      
  }
}

# draw plots
plot_kde206v207
plot_kde206v208

(a) KDEs based on ²⁰⁶Pb/²⁰⁴Pb vs ²⁰⁷Pb/²⁰⁴Pb.

(b) KDEs based on ²⁰⁶Pb/²⁰⁴Pb vs ²⁰⁸Pb/²⁰⁴Pb.

Figure 20.4: Binary plots with the KDEs of mining regions. Note that three different color gradients from light to dark represent respective intervals of 95%, 75%, and 50%.

20.4.5 Ternary diagram

This plotting method was first utilised by Cannon et al. (1961) who aimed to understand the principles of isotopic variations in ore lead. Raw isotopic ratios were expressed as relative abundances of ²⁰⁶Pb, ²⁰⁷Pb, ²⁰⁸Pb summed up to 100% by leaving out ²⁰⁴Pb. The transformed data were plotted as trilinear coordinates. The choice of represented masses was justified by the incapability of precisely determining the amount of ²⁰⁴Pb due to its low natural abundance of only 1.4% (uncertainties ~2.5%). These ternary diagrams were proposed as a solution to overcome the problem of analytically fundamentally biased data. However, they were rendered irrelevant with the advent of improved analytical techniques and the development of error correction models, which greatly increased precision (Taylor et al. 2015). A good example of using ternary diagrams is provided in (Hsu and Sabatini 2019).

The raw isotopic ratios are mathematically converted to three individual Pb compositions using the following equations:

\[ ^{206}Pb = \frac{\left(\frac{^{206}Pb}{^{204}Pb}\right) \cdot 100}{\left(\frac{^{206}Pb}{^{204}Pb}\right) + \left(\frac{^{207}Pb}{^{204}Pb}\right) + \left(\frac{^{208}Pb}{^{204}Pb}\right)} \]

\[ ^{207}Pb = \frac{\left(\frac{^{207}Pb}{^{204}Pb}\right) \cdot 100}{\left(\frac{^{206}Pb}{^{204}Pb}\right) + \left(\frac{^{207}Pb}{^{204}Pb}\right) + \left(\frac{^{208}Pb}{^{204}Pb}\right)} \]

\[ ^{208}Pb = \frac{\left(\frac{^{208}Pb}{^{204}Pb}\right) \cdot 100}{\left(\frac{^{206}Pb}{^{204}Pb}\right) + \left(\frac{^{207}Pb}{^{204}Pb}\right) + \left(\frac{^{208}Pb}{^{204}Pb}\right)} \]

Figure 20.5 displays a ternary scatter plot.

Show the code

## select region and columns with relative abundance of single isotopes 
data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}$", names(Roman_case), value = TRUE)))

# create plot 
plt <- plot_ly(data = data, a = ~Pb206, b = ~Pb207, c = ~Pb208,
               split = ~Region,  color = ~Region, symbol = ~Region, text = ~Region, 
               type = "scatterternary", 
               colors = colours, alpha = 0.75, symbols = symbols, 
               marker = list(line = list(color = 'rgb(0, 0, 0)', width = 1)),
               hovertemplate = paste('Region: %{text}<br>',
                                     '<sup>206</sup>Pb (%): %{a}<br>',
                                     '<sup>207</sup>Pb (%): %{b}<br>', 
                                     '<sup>208</sup>Pb (%): %{c}', 
                                     '<extra></extra>')
               ) |> 
  layout(ternary = list(
    sum = 100,
    aaxis = list(title = "<sup>206</sup>Pb (%)",
                 tickfont = list(size = 10),
                 tickangle = 60,
                 ticklen = 10,
                 min = min(data$Pb206)-0.05,
                 max = max(data$Pb206)+0.05
    ),
    baxis = list(title = "<sup>207</sup>Pb (%)",
                 tickfont = list(size = 10),
                 tickangle = -60,
                 ticklen = 10,
                 min = min(data$Pb207)-0.05,
                 max = max(data$Pb207)+0.05
    ),
    caxis = list(title = "<sup>208</sup>Pb (%)",
                 tickfont = list(size = 10),
                 ticklen = 10,
                 min = min(data$Pb208)-0.05,
                 max = max(data$Pb208)+0.05
                 )
    ), 
    legend = list(title = "", x = 0.75, y = 1.125)
  )

# draw plot
plt

Figure 20.5: A ternary plot of lead isotope data visualising the same dataset.

20.4.6 Ternary diagram with the kernel density estimation

KDE contour plots like in Figure 20.4 can also be generated for ternary diagrams (Figure 20.6). This helps us to better visualise the isotopic distribution of ore populations when many data are presented. However, the ternary KDE, in a sense, is not equal to the three-dimensional KDE. It is rather a regular KDE that is truncated to the ternary triangle.

Show the code

# create plot
plt <- ggtern(data = data, aes(x = Pb207, y = Pb206, z = Pb208, color = Region, fill = Region))+
  geom_density_tern(data = data[data$Region != "Roman object", ], n = 1000, bins=10, alpha = 1) +
  geom_point(data = data[data$Region == "Roman object", ], size = 2, alpha = 0.5) +
  scale_colour_manual(values = colours) +
  scale_T_continuous (name = expression({}^206*"Pb (%)"), limits = c(0.258,0.252)) +
  scale_L_continuous (name = expression({}^207*"Pb (%)"), limits = c(0.213,0.219)) +
  scale_R_continuous (name = expression({}^208*"Pb (%)"), limits = c(0.529,0.535)) +
  theme_bw() +
  theme(legend.position.inside = c(0.95,0.95),
        legend.justification = c(1,1),
        legend.box.just = 'left',
        legend.text = element_text(size = 12),
        legend.key = element_rect(fill = NA, colour=NA),
        plot.margin = unit(c(0,0,0,0), "mm")
        )

# draw plot
plt

Figure 20.6: Ternary plot with the KDEs of mining regions. Note that the contours in the KDEs represent different intervals.

20.4.7 Three-dimensional plot

The use of any single bivariate plot is insufficient for provenancing and is visually confusing when the ratios overlap. Therefore, additional diagrams are needed to show other combinations of isotopes. Three-dimensional plots represent the distribution of data in a three dimensional space (Figure 20.7) which has a higher discrimination power and is therefore better suited for provenance studies. The downside is that it is inherently difficult to read a 3D diagram and, therefore, a rotatable version is highly recommended.

Show the code

# extract data
data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}\\.Pb204", names(Roman_case), value = TRUE)))

# Create plot
plot_ly(data, 
        x = ~Pb206.Pb204, 
        y = ~Pb207.Pb204, 
        z = ~Pb208.Pb204, 
        color = ~Region, 
        colors = colours, 
        symbol = ~Region, 
        symbols = c("diamond", "triangle-up", "square", "cross"),
        text = ~Region,
        hovertemplate = paste('Region: %{text}<br>',
                              'x: %{x}<br>',
                              'y: %{y}<br>', 
                              'z: %{z}', 
                              '<extra></extra>')
        ) |>
  layout(scene = list(xaxis = list(title = '<sup>206</sup>Pb/<sup>204</sup>Pb'),
                      yaxis = list(title = '<sup>207</sup>Pb/<sup>204</sup>Pb'),
                      zaxis = list(title = '<sup>208</sup>Pb/<sup>204</sup>Pb'),
                      aspectmode='cube'
                      )
         )

Figure 20.7: The same dataset visualised as a three-dimensional plot.

20.4.8 Three-dimensional plot with kernel density estimation

This plotting method applies kernel density estimation to a three-dimensional diagram (Figure 20.8). It can help to delineate reference datasets with which targeted artefacts can be compared. Beardah and Baxter (1999) pioneered the application of a three-dimensional kernel plot in Pb isotope studies and suggested a sample size of 20 as an acceptable value. However, they also realised that larger sample sizes, ranging from 40 to 60, would be necessary if the population from which the sample is drawn is not normally distributed. To construct the 3D kernel plot using R, we modified the code from Ma et al. (2022).

Show the code

# select individual regions/assemblages
regions <- unique(data$Region)

# define levels, to be read as 100% - value, i.e. '5%' indicates 95% contour
levels <- c(95, 75, 50) 

# initiate data and plot
data_kde3d <- vector(mode = "list", length = length(regions)) 
names(data_kde3d) <- regions

par3d(windowRect = c(0, 0, 1200, 800))

# calculate KDEs for each region 
for (i in regions) {
  
  # subset data
  data_region <- data[data$Region == i, ]
  data_region$Region <- NULL # remove region for computation of KDEs
  
  # compute KDE
  data_region.pi <- Hpi(data_region, binned = TRUE) ## b is a matrix of x,y,z points
  data_kde3d[[i]] <- kde(data_region, H = data_region.pi)
  
  # plot KDE
  if (i != "Roman object") {
    plot(data_kde3d[[i]], display = "rgl", cont = levels, 
         col = colours[i], col.cont = colours[i], xlab = "", ylab = "", zlab = "", 
         add = TRUE, drawpoints = FALSE, box=FALSE, axes=FALSE
         )
  } else{
    plot(data_kde3d[[i]], display = "rgl", cont = 0, 
         col.pt = colours[i], size = 6, pch = 3, alpha = 1, xlab = "", ylab = "", zlab = "",
         add = TRUE, drawpoints = TRUE, box = TRUE, axes = FALSE
         )
  }
}

# decorate plot
title3d(xlab = expression({}^206*"Pb/"*{}^204*"Pb"), 
        ylab = expression({}^207*"Pb/"*{}^204*"Pb"), 
        zlab = expression({}^208*"Pb/"*{}^204*"Pb")
        )
bbox3d(color=c("grey","black"), emission="white", specular="lightgrey", shininess=5, alpha=0.8 )

# add legend
legend3d("topright", legend = regions, pch = 16, col = colours, cex = 2, inset = 0)

Figure 20.8: Three-dimensional plot with the KDEs of the mining regions.

20.5 Self check

Nowadays, which pairs of Pb isotopes are better suited to discriminate the isotopic ratios of artefacts and ore samples?
What statistical assumptions are appropriate to describe the distribution of Pb isotopic data in an assemblage?
What are the pros and cons when it comes to a three-dimensional diagram?

20.6 Further reading

Albarede F, Blichert-Toft J, Gentelli L, Milot J, Vaxevanopoulos M, Klein S, Westner KJ, Birch T, Davis G, Callataÿ F de (2020) A miner’s perspective on Pb isotope provenances in the Western and Central Mediterranean. J. Archaeol. Sci. 121:105194. https://doi.org/10.1016/j.jas.2020.105194
Blichert-Toft J, Delile H, Lee C-T, Stos-Gale Z, Billström K, Andersen T, Hannu H, Albarède F (2016) Large-scale tectonic cycles in Europe revealed by distinct Pb isotope provinces. Geochem. Geophys. Geosyst. 17:3854–3864. https://doi.org/10.1002/2016GC006524
Hsu Y-K, Sabatini BJ (2019) A geochemical characterization of lead ores in China: An isotope database for provenancing archaeological materials. PLoS ONE 14:e0215973. https://doi.org/10.1371/journal.pone.0215973

Albarède F, Blichert-Toft J, Gentelli L, et al (2020) A miner’s perspective on Pb isotope provenances in the Western and Central Mediterranean. Journal of Archaeological Science 121:105194. https://doi.org/10.1016/j.jas.2020.105194

Albarède F, Desaulty A-M, Blichert-Toft J (2012) A geological perspective on the use of Pb isotopes in Archaeometry. Archaeometry 54:853–867. https://doi.org/10.1111/j.1475-4754.2011.00653.x

Baxter MJ, Beardah CC, Wright RVS (1997) Some Archaeological Applications of Kernel Density Estimates. Journal of Archaeological Science 24:347–354. https://doi.org/10.1006/jasc.1996.0119

Baxter MJ, Gale NH (1998) Testing for multivariate normality via univariate tests: A case study using lead isotope ratio data. Journal of Applied Statistics 25:671–683. https://doi.org/10.1080/02664769822891

Beardah CC, Baxter MJ (1999) Three-Dimensional Data Display Using Kernel Density Estimates. In: Barceló JA, Briz I, Vila A (eds) New Techniques for Old Times. CAA98. Computer Applications and Quantitative Methods in Archaeology. Archaeopress, Oxford, pp 163–170

Bode M (2008) Archäometallurgische Untersuchungen zur Blei-/Silbergewinnung im Germanien der frühen Römischen Kaiserzeit. {PhD thesis}, Westfälische Wilhelms-Universität Münster

Cannon RS, Pierce AP, Antweiler JC, Buck KL (1961) The data of lead isotope geology related to problems of ore genesis. Economic Geology 56:1–38. https://doi.org/10.2113/gsecongeo.56.1.1

Ellam RM (2010) The graphical presentation of lead isotope data for environmental source apportionment. The Science of the total environment 408:3490–3492. https://doi.org/10.1016/j.scitotenv.2010.03.037

Hsu Y-K, Rawson J, Pollard AM, et al (2018) Application of Kernel Density Estimates to Lead Isotope Compositions of Bronzes from Ningxia, North-West China: Application of KDEs to Lead Isotope Compositions. Archaeometry 60:128–143. https://doi.org/10.1111/arcm.12347

Hsu Y-K, Sabatini BJ (2019) A geochemical characterization of lead ores in China: An isotope database for provenancing archaeological materials. PLoS ONE 14:e0215973. https://doi.org/10.1371/journal.pone.0215973

Ma Q, Pollard AM, Yu Y, et al (2022) Laser ablation inductively coupled plasma mass spectrometry analysis of potash and m-Na-Al glasses in China- using Kernel methods for trace element analysis. Heritage Science 10:253. https://doi.org/10.1186/s40494-022-00651-3

Scaife B, Budd P, McDonnell JG, Pollard AM (1999) Lead Isotope Analysis, Oxhide Ingots and the Presentation of Scientific Data in Archaeology. In: Young S, Pollard AM, Budd P, Ixer RA (eds) Metals in antiquity: Proceedings of the International symposium at Harvard University, 10 to 13 September 1997. Archaeopress, Oxford, pp 122–133

Stos-Gale ZA, Maliotis G, Gale NH, Annetts N (1997) Lead isotope characteristics of the Cyprus copper ore deposits applied to provenance studies of copper oxhide ingots. Archaeometry 39:83–123. https://doi.org/10.1111/j.1475-4754.1997.tb00792.x

Taylor RN, Ishizuka O, Michalik A, et al (2015) Evaluating the precision of Pb isotope measurement by mass spectrometry. Journal of Analytical Atomic Spectrometry 30:198–213. https://doi.org/10.1039/C4JA00279B

--- title: Visualisation of lead isotope data chapter-image: chapter-images/chap20.png author: - name: Yiu-Kang Hsu orcid: 0000-0002-2439-4863 corresponding: true email: Yiu-Kang.Hsu@bergbaumuseum.de - name: Thomas Rose orcid: 0000-0002-8186-3566 url: https://copper-smelting.com/ lightbox: true --- ## Learning objective By the end of this unit you are aware about the benefits and limitations of the different ways lead isotope data can be visualised. ## Prior knowledge This section assumes familiarity with the prior learning materials on Pb isotope geochemistry, esp. [Chapter 3](isotope_system.qmd). ## Material This section discusses different types of plots. Interactive examples of these plots allows you to explore their suitability for the different research questions on your own. ## Learning content Because lead isotopes are represented by three independent ratios, e.g., ^206^Pb/^204^Pb, ^207^Pb/^204^Pb, and ^208^Pb/^204^Pb, they can be visualised in a three dimensional geometric space. However, often only two dimensions can be represented in one plot. In addition, different ways of optically grouping and highlighting groups within the Pb isotope space exist. It is important to keep in mind that every plot has its pros and cons and there is no general consensus which is the best presentation. All plots were created in [R](https://www.r-project.org) with support of different packages. The source of the respective chunks is provided in collapsible sections. The entire source can be copied from the raw version of this chapter through the code tools next to the chapter heading. ```{r} #| label: load-pkgs #| code-summary: "Load R packages and datasets" #| message: false #| warning: false library(ggplot2) # creating static plots library(plotly) # creating interactive plots library(ggtern) # ternary plots library(ks) # calculation of kernel densities library(rgl) # creating interactive 3D plots library(viridisLite) # colour blind-friendly colour schemes setupKnitr(autoprint = TRUE) ``` The same subset of a real archaeological dataset will be used as example to demonstrate how this data can be treated in different ways. The dataset is available in [the GitHub repository of this book](https://github.com/archmetalDBM/GlobaLID-Edu/tree/main/example_dataset) and was originally published by [@Bode.2008]. To avoid overloading of the plots, the subset is restricted to the archaeological objects and reference data from the region "Northern Eifel", "Sauerland", and "Bergisches Land", all located in the West of Germany. ```{r} #| label: load-data #| code-summary: "Import and filter data" #| message: false Roman_case <- read.csv("example_dataset/Roman.csv", encoding="UTF-8", header = TRUE) |> subset(Region %in% c("Northern Eifel", "Sauerland", "Bergisches Land", "Roman object")) ``` In addition, the plots share several stylistic elements, such as the colour scheme. Therefore, they can be defined for all plots. ```{r} #| label: set-styling #| code-summary: "Define common style elements" #| message: false # define plotting order of regions Roman_case$Region <- factor(Roman_case$Region, levels = c('Northern Eifel', 'Sauerland', 'Bergisches Land', 'Roman object')) # define colours colours <- c(viridis(3), "black")|> setNames(c('Northern Eifel', 'Sauerland', 'Bergisches Land', 'Roman object')) # define symbols # in R, symbols are defined by numbers, see e.g., https://cran.r-project.org/web/packages/ggplot2/vignettes/ggplot2-specs.html#point symbols <- c(23,24,22,3) |> setNames(c('Northern Eifel', 'Sauerland', 'Bergisches Land', 'Roman object')) ``` ### Binary scatter plot The bi-plot (@fig-plotlybinary) is by far the most common option to display lead isotope data. Since there are four isotopes of Pb, twelve combinations of isotopic ratios can be derived. The use of paired ratios depends on the instruments used and the scientific disciplines of the studies. In the early days, Pb isotopic ratios were often reported based on ^206^Pb-based ratios as ^204^Pb could not be measured precisely. However, in the 2000s, the advent of the multi-collector mass spectrometer (MC-ICP-MS) and the double- or triple-spiked technique created a huge amount of Pb isotope data with precisely measured ^204^Pb. Conventionally, environmental science tends to use the ratios based on ^206^Pb, which however generates plots with linear patterns and thus a low discrimination power [@Ellam.2010]. In geological literature, ratios based on ^204^Pb are commonplace which enable a better visualisation of system closure time (or model age) and U-Th-Pb composition (or µ and κ) of parental source(s) [@Albarede.2012]. However, it has to be kept in mind that all two-dimensional plots incompletely represent a dataset. All twelve combination plots are suggested to be tested to view the full isotopic extent of ore deposits [@Albarede.2020]. Ideally, the Pb isotopic ratios should be considered in a three-dimensional space. ```{r} #| label: fig-plotlybinary #| fig-cap: "A binary plot of Roman lead artefacts in comparison with ore districts in Germany." #| message: false #| code-summary: "Show the code" # select the region and all columns with isotope ratios (= all columns starting with Pb followed by three characters and a period as replacement character for the slash) data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}\\.", names(Roman_case), value = TRUE))) colnames(data) <- c('Region', '204Pb/206Pb', '207Pb/206Pb', '208Pb/206Pb', '206Pb/204Pb', '207Pb/204Pb', '208Pb/204Pb') # create plot plt <- plot_ly(data = data, x = ~`206Pb/204Pb`, y = ~`207Pb/204Pb`, split = ~Region, color = ~Region, symbol = ~Region, text = ~Region, colors = colours, alpha = 0.75, symbols = symbols, marker = list(line = list(color = 'rgb(0, 0, 0)', width = 1)), hovertemplate = paste('Region: %{text} ', 'x: %{x} ', 'y: %{y}', '<extra></extra>') ) # set-up drop down menus for the x and y variables # function to create the drop down menu elements # adapted from a function provided on https://stackoverflow.com/questions/75596474/add-a-legend-to-a-plotly-plot-with-drop-down-menus-in-r btn <- function(xLoc, m, data, xOry) { # x btn position, method, names, x or y list(type = "list", x = xLoc, xanchor = "left", y = 1.2, buttons = lapply( names(data[[1]])[-1], function(k) { # iterate over names switch (xOry, x = {args <- list( # a list for each trace (each legend entry) list(x = list(data[[1]][, k], data[[1]][, k], data[[1]][, k], data[[1]][, k])), # 4 lists, 4 traces list(xaxis = list(title = k))) }, y = {args = list( # a list for each trace (each legend list(x = list(data[[1]][, k], data[[1]][, k], data[[1]][, k], data[[1]][, k])), # 4 lists, 4 traces list(yaxis = list(title = list(text = k))))}, stop("xOrY must be either 'x' or 'y'.") ) list(method = m, label = k, args = args) } ) ) } # split regions into separate datasets regions <- split(data, data$Region) # add dropdown menus to plot plt <- plt %>% layout( xaxis = list(title = "206Pb/204Pb"), yaxis = list(title = "207Pb/204Pb"), updatemenus = list(btn(.2, "update", regions, "x"), btn(.8, "update", regions, "y")), showlegend = TRUE, annotations = list( list( text = "X-Axis:", x=0.04, y=1.18, xref='paper', yref='paper',xanchor = "left", showarrow=FALSE ), list( text = "Y-Axis:", x=0.63, y=1.18, xref='paper', yref='paper',xanchor = "left", showarrow=FALSE ) ) ) # draw plot plt ``` ### Bi-plot using geological-informed parameters Instead of isotopic ratios, @Albarede.2012 advocate the use of calculated geological model parameters, namely the model age (T), U/Pb (μ), and Th/U (κ) to discriminate potential ore sources in provenance studies (@fig-plotlygeobinary). As shown in [chapter 3](isotope_system.qmd), ^206^Pb, ^207^Pb, and ^208^Pb are generated by radioactive decay of their parental isotopes ^238^U, ^235^U, and ^232^Th, respectively. We can therefore calculate the model age, ^238^U/^204^Pb and ^232^Th/^238^U from the Pb isotope ratios determined for a given sample using the equations provided in @Albarede.2012 or any other of the Pb isotope models mentioned in chapter 3 by using, e.g., an [R script](https://github.com/archmetalDBM/GlobaLID-database/blob/main/calculate_model_ages.R). ```{r} #| label: fig-plotlygeobinary #| fig-cap: "A binary plot of Roman lead artefacts using geological parameters of model age (T), U/Pb (μ), and Th/U (κ)." #| message: false #| code-summary: "Show the code" # select columns Age, Mu, Kappa, Region and rename them data <- subset(Roman_case, select=c('Region', 'Age', 'Mu', 'Kappa')) |> na.omit() colnames(data) <- c('Region', 'Age','U/Pb','Th/U') # create plot plt <- plot_ly(data = data, x = ~Age, y = ~`U/Pb`, split = ~Region, color = ~Region, symbol = ~Region, text = ~Region, colors = colours, alpha = 0.75, symbols = symbols, marker = list(line = list(color = 'rgb(0, 0, 0)', width = 1)), hovertemplate = paste('Region: %{text} ', 'x: %{x} ', 'y: %{y}', '<extra></extra>') ) # set-up drop down menus for the x and y variables. # The function 'btn', defined in the code chunk above, is reused. # split regions into separate datasets regions <- split(data, data$Region) # add dropdown menus to plot plt <- plt %>% layout( xaxis = list(title = "Age"), yaxis = list(title = "U/Pb"), updatemenus = list(btn(.2, "update", regions, "x"), btn(.8, "update", regions, "y")), showlegend = TRUE, annotations = list( list( text = "X-Axis:", x=0.04, y=1.18, xref='paper', yref='paper',xanchor = "left", showarrow=FALSE ), list( text = "Y-Axis:", x=0.63, y=1.18, xref='paper', yref='paper',xanchor = "left", showarrow=FALSE ) ) ) # draw plot plt ``` ### Bi-plot with 90% confidence ellipse The increasing amount of available Pb isotope analyses resulted in the issue of how to effectively present overlapping data. One way to circumvent this problem is the creation of so-called "ore fields" that represent the extent of Pb isotopic variation given that the data follow a normal distribution. Some researchers have resorted to the 90% confidence ellipse (@fig-ellipse) as the reference for ore sourcing [@StosGale.1997]. This technique was severely criticised by @Baxter.1998, who demonstrated the non-normality of Pb isotope data in a number of instances. The 90% confidence ellipse is no longer considered suitable for representing an ore Pb isotopic population. ```{r} #| label: fig-ellipse #| fig-cap: "Binary scatter plots with 90% confidence ellipses around data points." #| fig-subcap: #| - "Ellipses based on ^206^Pb/^204^Pb vs ^207^Pb/^204^Pb." #| - "Ellipses based on ^206^Pb/^204^Pb vs ^208^Pb/^204^Pb." #| fig-height: 8 #| fig-width: 10 #| layout-ncol: 2 #| message: false #| code-summary: "Show the code" # select the region and all 204Pb-normalised isotope ratios data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}\\.Pb204", names(Roman_case), value = TRUE))) colnames(data) <- c('Region', '206Pb/204Pb', '207Pb/204Pb', '208Pb/204Pb') # create plot for 207Pb/204Pb vs 206Pb/204Pb p_207 <- ggplot(data, aes(x = `206Pb/204Pb`, y = `207Pb/204Pb`, shape = Region, fill = Region)) + geom_point(stroke = 0.2,size = 2.5, color = "black", alpha = 0.75) + scale_shape_manual(values = symbols) + scale_fill_manual(values = colours) + stat_ellipse(level = 0.9, type = "t", linewidth = 0.1, show.legend = NA) + labs(x = expression(""^"206"*"Pb/"^"204"*"Pb"), y = expression(""^"207"*"Pb/"^"204"*"Pb")) + theme_bw() p_207 <- ggplotly(p_207, tooltip = c("x", "y", "fill")) %>% layout(xaxis = list(title = "206Pb/204Pb"), yaxis = list(title = "207Pb/204Pb"), legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175) ) # create plot for 208Pb/204Pb vs 206Pb/204Pb p_208 <- ggplot(data, aes(x = `206Pb/204Pb`, y = `208Pb/204Pb`, shape = Region, fill = Region)) + geom_point(stroke = 0.2,size = 2.5, color = "black", alpha = 0.75) + scale_shape_manual(values = symbols) + scale_fill_manual(values = colours) + stat_ellipse(level = 0.9, type = "t", linewidth = 0.1, show.legend = NA) + labs(x = expression(""^"206"*"Pb/"^"204"*"Pb"), y = expression(""^"208"*"Pb/"^"204"*"Pb")) + theme_bw() p_208 <- ggplotly(p_208, tooltip = c("x", "y", "fill")) %>% layout(xaxis = list(title = "206Pb/204Pb"), yaxis = list(title = "208Pb/204Pb"), legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175) ) # draw plots p_207 p_208 ``` ### Binary plot with the kernel density estimation Due to the non-normality of Pb isotope data, both @Baxter.1997 and @Scaife.1999 advocated the use of more robust kernel density estimates (KDEs) to display the isotopic extent of orefields (@fig-2dKDE). KDEs are a non-parametric method to transform continuous data into a smoothed probability density function. KDEs offer three main advantages: 1. They do not assume the normality of data; 2. They can produce smoother distributions than conventional histograms, whose appearance is significantly affected by the choices of bin width and the start/end points of bins; and 3. They can represent data in a multidimensional space and enable users to effectively compare different datasets either graphically or mathematically. Given these advantages, the KDE method has become popular in recent publications of archaeological sciences [@Hsu.2018]. ```{r} #| label: fig-2dKDE #| fig-cap: "Binary plots with the KDEs of mining regions. Note that three different color gradients from light to dark represent respective intervals of 95%, 75%, and 50%." #| fig-subcap: #| - "KDEs based on ^206^Pb/^204^Pb vs ^207^Pb/^204^Pb." #| - "KDEs based on ^206^Pb/^204^Pb vs ^208^Pb/^204^Pb." #| fig-height: 8 #| fig-width: 10 #| message: false #| warning: false #| code-summary: "Show the code" #| layout-ncol: 2 # select individual regions/assemblages regions <- c("Northern Eifel", "Sauerland", "Bergisches Land") # define levels, to be read as 100% - value, i.e. '5%' indicates 95% contour levels <- c("5%", "25%", "50%") # initiate data data_kde206v207 <- vector(mode = "list", length = length(regions)) names(data_kde206v207) <- regions data_kde206v208 <- vector(mode = "list", length = length(regions)) names(data_kde206v208) <- regions for (i in regions) { # subset data data_region <- Roman_case[Roman_case$Region == i, ] # compute KDE data_region_kde206v207 <- ks::kde(data_region[c("Pb206.Pb204", "Pb207.Pb204")], compute.cont = TRUE,gridsize = 1024) data_region_kde206v208 <- ks::kde(data_region[c("Pb206.Pb204", "Pb208.Pb204")], compute.cont = TRUE,gridsize = 1024) # compute contour lines data_region_kde206v207_levels <- data_region_kde206v207$cont[levels] data_region_kde206v207_contour <- with(data_region_kde206v207, contourLines(x = eval.points[[1]], y = eval.points[[2]], z = estimate, levels = data_region_kde206v207_levels)) names(data_region_kde206v207_contour) <- paste0(i, seq_along(data_region_kde206v207_contour)) data_region_kde206v207_contour <- Map(c, data_region_kde206v207_contour, group = names(data_region_kde206v207_contour)) data_region_kde206v207_contour <- do.call("rbind", lapply(data_region_kde206v207_contour, data.frame)) data_region_kde206v207_contour$level_percent <- names(data_region_kde206v207_levels)[match(data_region_kde206v207_contour$level, data_region_kde206v207_levels)] data_region_kde206v208_levels <- data_region_kde206v208$cont[levels] data_region_kde206v208_contour <- with(data_region_kde206v208, contourLines(x = eval.points[[1]], y = eval.points[[2]], z = estimate, levels = data_region_kde206v208_levels)) names(data_region_kde206v208_contour) <- paste0(i, seq_along(data_region_kde206v208_contour)) data_region_kde206v208_contour <- Map(c, data_region_kde206v208_contour, group = names(data_region_kde206v208_contour)) data_region_kde206v208_contour <- do.call("rbind", lapply(data_region_kde206v208_contour, data.frame)) data_region_kde206v208_contour$level_percent <- names(data_region_kde206v208_levels)[match(data_region_kde206v208_contour$level, data_region_kde206v208_levels)] data_kde206v207[[i]] <- data_region_kde206v207_contour data_kde206v208[[i]] <- data_region_kde206v208_contour } # crunch list into data frame data_kde206v207 <- Map(cbind, data_kde206v207, Region = names(data_kde206v207)) data_kde206v207 <- do.call("rbind", data_kde206v207) data_kde206v208 <- Map(cbind, data_kde206v208, Region = names(data_kde206v208)) data_kde206v208 <- do.call("rbind", data_kde206v208) # plots plot_kde206v207 <- ggplot() + geom_point(data = Roman_case[Roman_case$Region == "Roman object", ], aes(x = Pb206.Pb204, y = Pb207.Pb204, colour = Region), shape = 3, size = 2.3) + geom_polygon(data = data_kde206v207, mapping = aes(x = x, y= y, group = group, colour = Region, fill = Region, alpha = level_percent)) + scale_colour_manual(name = NULL, limits = names(colours), values = colours) + scale_fill_manual(name = NULL, limits = names(colours), values = colours) + scale_alpha_discrete(range = c(0.1, 0.3), guide = "none") + scale_x_continuous(name = expression(""^"206"*"Pb/"^"204"*"Pb"), limits = c(18, 18.6), breaks = seq(18, 18.6, by = 0.1)) + scale_y_continuous(name = expression(""^"207"*"Pb/"^"204"*"Pb"), limits = c(15.5, 15.8), breaks = seq(15.5, 15.8, by = 0.05)) + guides(colour = guide_legend(position = "inside")) + theme_bw() + theme(panel.grid.minor = element_blank(), axis.text = element_text(size = 10), axis.title = element_text(size = 14), legend.position.inside = c(0.76, 0.11) ) plot_kde206v208 <- ggplot() + geom_point(data = Roman_case[Roman_case$Region == "Roman object", ], aes(x = Pb206.Pb204, y = Pb208.Pb204, colour = Region), shape = 3, size = 2.3) + geom_polygon(data = data_kde206v208, mapping = aes(x = x, y = y, group = group, colour = Region, fill = Region, alpha = level_percent)) + geom_point(data = Roman_case[Roman_case$Region == "Roman object", ], aes(x = Pb206.Pb204, y = Pb208.Pb204, colour = Region), shape = 3, size = 2.3) + scale_colour_manual(name = NULL, limits = names(colours), values = colours) + scale_fill_manual(name = NULL, limits = names(colours), values = colours) + scale_alpha_discrete(range = c(0.1, 0.3), guide = "none") + scale_x_continuous(name = expression(""^"206"*"Pb/"^"204"*"Pb"), limits = c(18, 18.6), breaks = seq(18, 18.6, by = 0.1)) + scale_y_continuous(name = expression(""^"208"*"Pb/"^"204"*"Pb"), limits = c(37.9, 38.8), breaks = seq(38, 39, by = 0.1)) + guides(colour = guide_legend(position = "inside")) + theme_bw() + theme(panel.grid.minor = element_blank(), axis.text = element_text(size = 10), axis.title = element_text(size = 14), legend.position.inside = c(0.76, 0.11) ) # convert to plotly plot_kde206v207 <- ggplotly(plot_kde206v207) %>% layout(xaxis = list(title = "206Pb/204Pb"), yaxis = list(title = "207Pb/204Pb"), legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175) ) plot_kde206v208 <- ggplotly(plot_kde206v208) %>% layout(xaxis = list(title = "206Pb/204Pb"), yaxis = list(title = "207Pb/204Pb"), legend = list(title = "", font = list(size = 10), orientation = "h", xanchor = "center", x = 0.5, y = 1.175) ) # modify legends for (i in seq_along(plot_kde206v207$x$data)) { if (!grepl("25%|Roman", plot_kde206v207$x$data[[i]]$name)) {plot_kde206v207$x$data[[i]]$showlegend <- FALSE} for (k in unique(Roman_case$Region)) { if (grepl(k, plot_kde206v207$x$data[[i]]$name)) { plot_kde206v207$x$data[[i]]$name <- k plot_kde206v207$x$data[[i]]$legendgroup <- k } } } for (i in seq_along(plot_kde206v208$x$data)) { if (!grepl("25%|Roman", plot_kde206v208$x$data[[i]]$name)) {plot_kde206v208$x$data[[i]]$showlegend <- FALSE} for (k in unique(Roman_case$Region)) { if (grepl(k, plot_kde206v208$x$data[[i]]$name)) { plot_kde206v208$x$data[[i]]$name <- k plot_kde206v208$x$data[[i]]$legendgroup <- k } } } # draw plots plot_kde206v207 plot_kde206v208 ``` ### Ternary diagram This plotting method was first utilised by @Cannon.1961 who aimed to understand the principles of isotopic variations in ore lead. Raw isotopic ratios were expressed as relative abundances of ^206^Pb, ^207^Pb, ^208^Pb summed up to 100% by leaving out ^204^Pb. The transformed data were plotted as trilinear coordinates. The choice of represented masses was justified by the incapability of precisely determining the amount of ^204^Pb due to its low natural abundance of only 1.4% (uncertainties \~2.5%). These ternary diagrams were proposed as a solution to overcome the problem of analytically fundamentally biased data. However, they were rendered irrelevant with the advent of improved analytical techniques and the development of error correction models, which greatly increased precision [@Taylor.2015]. A good example of using ternary diagrams is provided in [@Hsu.2019]. The raw isotopic ratios are mathematically converted to three individual Pb compositions using the following equations: $$ ^{206}Pb = \frac{\left(\frac{^{206}Pb}{^{204}Pb}\right) \cdot 100}{\left(\frac{^{206}Pb}{^{204}Pb}\right) + \left(\frac{^{207}Pb}{^{204}Pb}\right) + \left(\frac{^{208}Pb}{^{204}Pb}\right)} $$ $$ ^{207}Pb = \frac{\left(\frac{^{207}Pb}{^{204}Pb}\right) \cdot 100}{\left(\frac{^{206}Pb}{^{204}Pb}\right) + \left(\frac{^{207}Pb}{^{204}Pb}\right) + \left(\frac{^{208}Pb}{^{204}Pb}\right)} $$ $$ ^{208}Pb = \frac{\left(\frac{^{208}Pb}{^{204}Pb}\right) \cdot 100}{\left(\frac{^{206}Pb}{^{204}Pb}\right) + \left(\frac{^{207}Pb}{^{204}Pb}\right) + \left(\frac{^{208}Pb}{^{204}Pb}\right)} $$ @fig-ternaryscatter displays a ternary scatter plot. ```{r} #| label: fig-ternaryscatter #| fig-cap: "A ternary plot of lead isotope data visualising the same dataset." #| message: false #| code-summary: "Show the code" ## select region and columns with relative abundance of single isotopes data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}$", names(Roman_case), value = TRUE))) # create plot plt <- plot_ly(data = data, a = ~Pb206, b = ~Pb207, c = ~Pb208, split = ~Region, color = ~Region, symbol = ~Region, text = ~Region, type = "scatterternary", colors = colours, alpha = 0.75, symbols = symbols, marker = list(line = list(color = 'rgb(0, 0, 0)', width = 1)), hovertemplate = paste('Region: %{text} ', '206Pb (%): %{a} ', '207Pb (%): %{b} ', '208Pb (%): %{c}', '<extra></extra>') ) |> layout(ternary = list( sum = 100, aaxis = list(title = "206Pb (%)", tickfont = list(size = 10), tickangle = 60, ticklen = 10, min = min(data$Pb206)-0.05, max = max(data$Pb206)+0.05 ), baxis = list(title = "207Pb (%)", tickfont = list(size = 10), tickangle = -60, ticklen = 10, min = min(data$Pb207)-0.05, max = max(data$Pb207)+0.05 ), caxis = list(title = "208Pb (%)", tickfont = list(size = 10), ticklen = 10, min = min(data$Pb208)-0.05, max = max(data$Pb208)+0.05 ) ), legend = list(title = "", x = 0.75, y = 1.125) ) # draw plot plt ``` ### Ternary diagram with the kernel density estimation KDE contour plots like in @fig-2dKDE can also be generated for ternary diagrams (@fig-ternaryKDE). This helps us to better visualise the isotopic distribution of ore populations when many data are presented. However, the ternary KDE, in a sense, is not equal to the three-dimensional KDE. It is rather a regular KDE that is truncated to the ternary triangle. ```{r} #| label: fig-ternaryKDE #| fig-cap: "Ternary plot with the KDEs of mining regions. Note that the contours in the KDEs represent different intervals. " #| message: false #| warning: false #| code-summary: "Show the code" # create plot plt <- ggtern(data = data, aes(x = Pb207, y = Pb206, z = Pb208, color = Region, fill = Region))+ geom_density_tern(data = data[data$Region != "Roman object", ], n = 1000, bins=10, alpha = 1) + geom_point(data = data[data$Region == "Roman object", ], size = 2, alpha = 0.5) + scale_colour_manual(values = colours) + scale_T_continuous (name = expression({}^206*"Pb (%)"), limits = c(0.258,0.252)) + scale_L_continuous (name = expression({}^207*"Pb (%)"), limits = c(0.213,0.219)) + scale_R_continuous (name = expression({}^208*"Pb (%)"), limits = c(0.529,0.535)) + theme_bw() + theme(legend.position.inside = c(0.95,0.95), legend.justification = c(1,1), legend.box.just = 'left', legend.text = element_text(size = 12), legend.key = element_rect(fill = NA, colour=NA), plot.margin = unit(c(0,0,0,0), "mm") ) # draw plot plt ``` ### Three-dimensional plot The use of any single bivariate plot is insufficient for provenancing and is visually confusing when the ratios overlap. Therefore, additional diagrams are needed to show other combinations of isotopes. Three-dimensional plots represent the distribution of data in a three dimensional space (@fig-3dplot) which has a higher discrimination power and is therefore better suited for provenance studies. The downside is that it is inherently difficult to read a 3D diagram and, therefore, a rotatable version is highly recommended. ```{r} #| label: fig-3dplot #| fig-cap: "The same dataset visualised as a three-dimensional plot. " #| message: false #| code-summary: "Show the code" # extract data data <- subset(Roman_case, select = c('Region', grep("^Pb.{3}\\.Pb204", names(Roman_case), value = TRUE))) # Create plot plot_ly(data, x = ~Pb206.Pb204, y = ~Pb207.Pb204, z = ~Pb208.Pb204, color = ~Region, colors = colours, symbol = ~Region, symbols = c("diamond", "triangle-up", "square", "cross"), text = ~Region, hovertemplate = paste('Region: %{text} ', 'x: %{x} ', 'y: %{y} ', 'z: %{z}', '<extra></extra>') ) |> layout(scene = list(xaxis = list(title = '206Pb/204Pb'), yaxis = list(title = '207Pb/204Pb'), zaxis = list(title = '208Pb/204Pb'), aspectmode='cube' ) ) ``` ### Three-dimensional plot with kernel density estimation This plotting method applies kernel density estimation to a three-dimensional diagram (@fig-3dKDE). It can help to delineate reference datasets with which targeted artefacts can be compared. @Beardah.1999 pioneered the application of a three-dimensional kernel plot in Pb isotope studies and suggested a sample size of 20 as an acceptable value. However, they also realised that larger sample sizes, ranging from 40 to 60, would be necessary if the population from which the sample is drawn is not normally distributed. To construct the 3D kernel plot using R, we modified the code from @Ma.2022. ```{r} #| label: fig-3dKDE #| message: false #| warning: false #| code-summary: Show the code #| fig-cap: Three-dimensional plot with the KDEs of the mining regions. # select individual regions/assemblages regions <- unique(data$Region) # define levels, to be read as 100% - value, i.e. '5%' indicates 95% contour levels <- c(95, 75, 50) # initiate data and plot data_kde3d <- vector(mode = "list", length = length(regions)) names(data_kde3d) <- regions par3d(windowRect = c(0, 0, 1200, 800)) # calculate KDEs for each region for (i in regions) { # subset data data_region <- data[data$Region == i, ] data_region$Region <- NULL # remove region for computation of KDEs # compute KDE data_region.pi <- Hpi(data_region, binned = TRUE) ## b is a matrix of x,y,z points data_kde3d[[i]] <- kde(data_region, H = data_region.pi) # plot KDE if (i != "Roman object") { plot(data_kde3d[[i]], display = "rgl", cont = levels, col = colours[i], col.cont = colours[i], xlab = "", ylab = "", zlab = "", add = TRUE, drawpoints = FALSE, box=FALSE, axes=FALSE ) } else{ plot(data_kde3d[[i]], display = "rgl", cont = 0, col.pt = colours[i], size = 6, pch = 3, alpha = 1, xlab = "", ylab = "", zlab = "", add = TRUE, drawpoints = TRUE, box = TRUE, axes = FALSE ) } } # decorate plot title3d(xlab = expression({}^206*"Pb/"*{}^204*"Pb"), ylab = expression({}^207*"Pb/"*{}^204*"Pb"), zlab = expression({}^208*"Pb/"*{}^204*"Pb") ) bbox3d(color=c("grey","black"), emission="white", specular="lightgrey", shininess=5, alpha=0.8 ) # add legend legend3d("topright", legend = regions, pch = 16, col = colours, cex = 2, inset = 0) ``` ## Self check - Nowadays, which pairs of Pb isotopes are better suited to discriminate the isotopic ratios of artefacts and ore samples? - What statistical assumptions are appropriate to describe the distribution of Pb isotopic data in an assemblage? - What are the pros and cons when it comes to a three-dimensional diagram? ## Further reading - Albarede F, Blichert-Toft J, Gentelli L, Milot J, Vaxevanopoulos M, Klein S, Westner KJ, Birch T, Davis G, Callataÿ F de (2020) A miner's perspective on Pb isotope provenances in the Western and Central Mediterranean. J. Archaeol. Sci. 121:105194. <https://doi.org/10.1016/j.jas.2020.105194> - Blichert-Toft J, Delile H, Lee C-T, Stos-Gale Z, Billström K, Andersen T, Hannu H, Albarède F (2016) Large-scale tectonic cycles in Europe revealed by distinct Pb isotope provinces. Geochem. Geophys. Geosyst. 17:3854–3864. <https://doi.org/10.1002/2016GC006524> - Hsu Y-K, Sabatini BJ (2019) A geochemical characterization of lead ores in China: An isotope database for provenancing archaeological materials. PLoS ONE 14:e0215973. <https://doi.org/10.1371/journal.pone.0215973>