---
title: "pXRF Data Analysis and Visualization"
author: "Martin Valent"
date: "`r Sys.Date()`"
format:
  html:
    toc: true
    code-fold: true
execute:
  echo: true
---

# Introduction

This Quarto document contains R code for analyzing and visualizing pXRF data. The analysis includes copper and lead categorization, heatmaps, boxplots, and PCA analysis of elemental concentrations in drilled samples.

------------------------------------------------------------------------

# Required Libraries

```{r}
install.packages("pheatmap")
install.packages("dplyr")
install.packages("here")
install.packages("reshape2")
install.packages("ggplot2")
install.packages("compositions")
install.packages("ggbiplot")

library(pheatmap)
library(dplyr)
library(here)
library(reshape2)
library(ggplot2)
library(compositions)
library(ggbiplot)

```

##Copper categories pheatmap

```{r}
#Copper categoris - "Cate" dataset; for this code you need datasets "basic_pXRF" and "bronzy_pxrf" which are in RAR. file
basic_pXRF <- read.csv(here("basic_pXRF.csv"), sep = ",")
bronzy_pxrf <- read.csv(here("bronzy_pxrf.csv"), sep = ",")

# Work on a copy of original dataset
basic_pXRF_cate <- basic_pXRF  #Cate > copper categories of surface pXRF measurements; 91 objects

# Define copper categories
copper_categories <- data.frame(
  Cu_category = 1:16,
  As_present = c(0,1,0,0,  0,1,0,0,  1,0,1,1,  0,1,1,1),
  Sb_present = c(0,0,1,0,  0,1,1,0,  0,1,0,1,  1,1,0,1),
  Ag_present = c(0,0,0,1,  0,0,1,1,  1,0,0,1,  1,0,1,1),
  Ni_present = c(0,0,0,0,  1,0,0,1,  0,1,1,0,  1,1,1,1)
)

#Remove unwanted elements
basic_pXRF_cate <- basic_pXRF_cate[, !colnames(basic_pXRF_cate) %in% c("Ti", "Fe", "Sn")]

# Assign copper categories and update 'urcenie' - typology
basic_pXRF_cate$urcenie <- bronzy_pxrf$urcenie
if (all(rownames(basic_pXRF_cate) %in% bronzy_pxrf$ID_nalezu)) {
  print("All row names match ID_nalezu values.")
} else {
  print("There are mismatches between row names and ID_nalezu.")
}

# Remove iron slag
basic_pXRF_cate <- basic_pXRF_cate[!basic_pXRF_cate$ID %in% c("28", "132"), ]

# Apply thresholding on the original data to generate binary variables for each element
threshold <- 0.1
basic_pXRF_cate$As_present <- ifelse(basic_pXRF_cate$As > threshold, 1, 0)
basic_pXRF_cate$Sb_present <- ifelse(basic_pXRF_cate$Sb > threshold, 1, 0)
basic_pXRF_cate$Ag_present <- ifelse(basic_pXRF_cate$Ag > threshold, 1, 0)
basic_pXRF_cate$Ni_present <- ifelse(basic_pXRF_cate$Ni > threshold, 1, 0)

# Merge artefact data with copper category definitions
basic_pXRF_cate <- merge(basic_pXRF_cate, copper_categories, 
                         by = c("As_present", "Sb_present", "Ag_present", "Ni_present"), 
                         all.x = TRUE)

# Prepare data for heatmap
basic_pXRF_cate$urcenie[basic_pXRF_cate$urcenie == "nástroje"] <- "Tools"
basic_pXRF_cate$urcenie[basic_pXRF_cate$urcenie == "neurčené"] <- "Undefined"
basic_pXRF_cate$urcenie[basic_pXRF_cate$urcenie == "čepele"] <- "Blades"
basic_pXRF_cate$urcenie[basic_pXRF_cate$urcenie == "zvyšky bronzovej metalurgickej produkcie"] <- "Bronzeworking waste"
basic_pXRF_cate$urcenie[basic_pXRF_cate$urcenie == "zvyšky železnej metalurgickej produkcie"] <- "Iron Slag"
basic_pXRF_cate$urcenie[basic_pXRF_cate$urcenie == "šperky a súčasti kroja"] <- "Jewellery"
basic_pXRF_cate$urcenie[basic_pXRF_cate$urcenie == "militáriá"] <- "Weapons"

copper_categories_names <- c("1" = "none", "2" = "As", "3" = "Sb", "4" = "Ag", "6" = "As+Sb", "9" = "As+Ag", "10" = "Sb+Ni", "11" = "As+Ni", "12" = "As+Sb+Ag", "14" = "As+Sb+Ni", "16" = "As+Sb+Ag+Ni")
heatmap_cate <- table(basic_pXRF_cate$Cu_category, basic_pXRF_cate$urcenie)
rownames(heatmap_cate) <- copper_categories_names[as.character(rownames(heatmap_cate))]
heatmap_cate

# Generate heatmap for copper categories
pheatmap(heatmap_cate,
         display_numbers = TRUE,
         number_format = "%.f",
         number_fontface = "bold",
         main = "Heatmap of Copper Categories",
         scale = "none",
         clustering_distance_rows = "euclidean",
         clustering_distance_cols = "euclidean",
         color = colorRampPalette(c("white", "yellow"))(100),
         legend = FALSE,
         angle = 45,
         cluster_rows = FALSE,
         cluster_cols = FALSE)
```

##Lead categories pheatmap

```{r}
#Lead frequency - "Pablo" dataset; for this code you need datasets "basic_pXRF" and "bronzy_pxrf" which are in RAR. file

# Work on a copy of original dataset
basic_pxrf_pablo <- basic_pXRF #Pablo > lead volume in surface pXRF measurements; 91 objects

#Define lead categories
Pb_category <- case_when(
  basic_pxrf_pablo$Pb < 0.1  ~ 0,
  basic_pxrf_pablo$Pb >= 0.1 & basic_pxrf_pablo$Pb < 1  ~ 1,
  basic_pxrf_pablo$Pb >= 1 & basic_pxrf_pablo$Pb < 5  ~ 2,
  basic_pxrf_pablo$Pb >= 5 & basic_pxrf_pablo$Pb < 10  ~ 3,
  basic_pxrf_pablo$Pb >= 10  ~ 4
)

# Exclude unwanted elements (Ti, Fe, Sn)
basic_pxrf_pablo <- basic_pxrf_pablo[, !colnames(basic_pxrf_pablo) %in% c("Ti", "Fe", "Sn")]


# Assign Pb categories and update 'urcenie' - typology
annotation_pb <- data.frame(Pb_category = as.factor(Pb_category))
rownames(annotation_pb) <- basic_pxrf_pablo$ID
annotation_pb$Pb_category <- as.character(annotation_pb$Pb_category)
basic_pxrf_pablo$Pb_category <- annotation_pb$Pb_category
basic_pxrf_pablo$urcenie <- bronzy_pxrf$urcenie

if (all(rownames(basic_pxrf_pablo) %in% bronzy_pxrf$ID_nalezu)) {
  print("All row names match ID_nalezu values.")
} else {
  print("There are mismatches between row names and ID_nalezu.")
}

#Remove discovered iron slag remains
basic_pxrf_pablo <- basic_pxrf_pablo[!basic_pxrf_pablo$ID %in% c("28", "132"), ]

#Addition of typological classification
basic_pxrf_pablo$urcenie[basic_pxrf_pablo$urcenie == "nástroje"] <- "Tools"
basic_pxrf_pablo$urcenie[basic_pxrf_pablo$urcenie == "neurčené"] <- "Undefined"
basic_pxrf_pablo$urcenie[basic_pxrf_pablo$urcenie == "čepele"] <- "Blades"
basic_pxrf_pablo$urcenie[basic_pxrf_pablo$urcenie == "zvyšky bronzovej metalurgickej produkcie"] <- "Bronzeworking waste"
basic_pxrf_pablo$urcenie[basic_pxrf_pablo$urcenie == "zvyšky železnej metalurgickej produkcie"] <- "Iron Slag"
basic_pxrf_pablo$urcenie[basic_pxrf_pablo$urcenie == "šperky a súčasti kroja"] <- "Jewellery"
basic_pxrf_pablo$urcenie[basic_pxrf_pablo$urcenie == "militáriá"] <- "Weapons"

# Prepare data for lead heatmap
heatmap_pablo <- table(basic_pxrf_pablo$Pb_category, basic_pxrf_pablo$urcenie)
pb_category_names <- c("0" = "Pb = <0,1%", "1" = "Pb = 0,1-1%", "2" = "Pb = 1-5%", "3" = "Pb = 5-10%", "4" = "Pb = >10%")
rownames(heatmap_pablo) <- pb_category_names[as.character(rownames(heatmap_pablo))]
heatmap_pablo

# Generate lead heatmap
pheatmap(heatmap_pablo,
         display_numbers = TRUE,
         number_format = "%.f",
         number_fontface = "bold",
         fontsize_number = 12,
         main = "Heatmap of Lead values",
         scale = "none",
         clustering_distance_rows = "euclidean",
         clustering_distance_cols = "euclidean",
         color = colorRampPalette(c("white", "grey"))(100),
         legend = FALSE,
         angle = 45,
         cluster_rows = FALSE,
         cluster_cols = FALSE)
```

##Boxplot with trace elements values in drilled samples

```{r}
# Ensure the 'reshape2' package is loaded
library(reshape2)

#This code requires "basic_xrf" dataset
basic_xrf <- read.csv(here("basic_xrf.csv"), sep = ",")

# Reshape the wide data (drilled samples) to long format for the boxplot
basic_xrf_box_prep <- basic_xrf

# Remove elements that are not required for the boxplot
basic_xrf_box_prep$Cu <- NULL
basic_xrf_box_prep$Co <- NULL
basic_xrf_box_prep$Zn <- NULL
basic_xrf_box_prep$Fe <- NULL
basic_xrf_box_prep$Au <- NULL

# Convert the dataset from wide to long format
basic_xrf_long <- melt(basic_xrf_box_prep, 
                       id.vars = "ID", 
                       variable.name = "Element", 
                       value.name = "Value")

# Create the boxplot for the measured values of the remaining elements
boxplot(Value ~ Element, data = basic_xrf_long,
        xlab = "", ylab = "Measured value",
        main = "",
        col = "lightblue", border = "darkblue")
```

##PCA analysis on drilled samples

```{r}
#Code for PCA analysis on drilled samples--------

# Load required libraries
library("here")
library("ggplot2")
library("compositions")
library("ggbiplot")

# Ensure basic_xrf has proper rownames by writing and re-reading
basic_xrf <- read.csv(here("basic_xrf.csv"))
rownames(basic_xrf) <- basic_xrf$ID
basic_xrf$ID <- NULL
basic_xrf_merged2 <- read.csv(here("basic_xrf_merged2.csv"), sep = ",")
numerical_xrf_data <- basic_xrf[, sapply(basic_xrf, is.numeric)] 
numerical_xrf_data <- numerical_xrf_data[!rownames(numerical_xrf_data)%in% c("28", "132"), ] 
numerical_xrf_data$Au <- NULL
numerical_xrf_data$Co <- NULL
numerical_xrf_data$Fe <- NULL
small_constant <- min(numerical_xrf_data[numerical_xrf_data > 0])*0.5
numerical_xrf_data[numerical_xrf_data == 0] <- small_constant

row_totals <- rowSums(numerical_xrf_data)
pseudo_total <- 100
rescaled_xrf_data <- sweep(numerical_xrf_data, 1, row_totals/pseudo_total, FUN = "/")
row_totals <- rowSums(rescaled_xrf_data)
clr_xrf_data <- clr(rescaled_xrf_data)

clr_xrf_data_matrix <- as.matrix(clr_xrf_data)
clr_xrf_df <- as.data.frame(clr_xrf_data_matrix)

pca_normalised <- prcomp(clr_xrf_df, center = TRUE, scale. = TRUE)

#Eigenvalues for CLR transformation
eigenvalues <- pca_normalised$sdev^2
eigenvalues_clr <- round(eigenvalues, 4)
eigenvalues_clr <- format(eigenvalues, scientific = FALSE)

#Loadings for CLR transformation
loadings_clr <- pca_normalised$rotation

#visualization
pca_scores <- pca_normalised$x

pca_scores <- as.data.frame(pca_scores)
pca_scores$PC1 <- as.numeric(as.character(pca_scores$PC1))
pca_scores$PC2 <- as.numeric(as.character(pca_scores$PC2))
pca_scores$PC3 <- as.numeric(as.character(pca_scores$PC3))

urcenie_xrf <- basic_xrf_merged2 #dataframe for "Urcenie" column addition
urcenie_xrf$ID <- as.character(urcenie_xrf$ID)
urcenie_xrf <- urcenie_xrf[!urcenie_xrf$ID %in% c(28, 132), ]
urcenie_xrf$urcenie <- gsub("Bronze Slag", "Bronzeworking waste", urcenie_xrf$urcenie)


if (identical(rownames(pca_scores), as.character(urcenie_xrf$ID))) {
  pca_scores <- cbind(pca_scores, urcenie = urcenie_xrf$urcenie)
  message("Binding is correct: Row names match 'ID' column values.")
} else {
  message("Warning: Row names in pca_scores do not match 'ID' column values.")
}

eigenvalues_clr
loadings_clr

ggbiplot(pca_normalised, obs.scale = 1, var.scale = 1, choices = c(1, 2), alpha = 0) +
  geom_point(data = pca_scores,
             aes(x = PC1, y = PC2, color = urcenie),
             size = 3) +
  geom_text(aes(label = rownames(pca_scores)), vjust = -1, hjust = 0.5, color = "black", size = 3)+
  scale_color_manual(values = rainbow(length(unique(pca_scores$urcenie)))) +
  labs(color = "Type",
       title = "PCA Biplot - drilled samples")+
  theme_minimal()
```