---
title: "Pilotstudie Farchant - Datenkontrolle, -aufbereitung und -auswertung"
subtitle: "Basierend auf den Daten des Bruker Tracer 5g No. 900F3980"
author: 
  name: Michaela Schauer
  orcid: 0000-0001-7514-7187
  affiliations:
    - Department für Kulturwissenschaften, Geschwister-Scholl-Platz 1, Ludwig-Maximilians-Universität München, München 80539, Deutschland
    - Vienna Institute for Archaeological Science (VIAS), Universität Wien, Franz-Klein-Gasse 1, 1190 Wien, Österreich
    - Research Network Human Evolution and Archaeological Sciences (HEAS), Universität Wien, Djerassiplatz 1, 1030 Wien, Österreich
date: "`r Sys.Date()`" 

toc: true
toc-depth: 5
toc-title: Contents
number-sections: true
callout-icon: false

knitr:
  opts_chunk: 
    warning: false
    message: false
    
editor: visual
license: CC BY

bibliography: R-Zitate_Farchant.bib

format:
  html:
    toc-location: left
    embed-resources: true
    code-overflow: wrap
---

\newpage

## Grundsätzliches

Dieses Quarto Skript (R Quarto v.1.5.55) [@Allaire.2024] wurde mit R v. 4.4.0 [@RCoreTeam.2024] und RStudio v. 2024.04.2 [@RStudioTeam.2024] erstellt. Es werden folgende R-Packages verwendet: basictabler [@Bailiss.2021] - cowplot [@Wilke.2024] - data.table [@Barrett.2024] - dplyr [@Wickham.2023c] - forcats [@Wickham.2023] - ggplot 2 [@Wickham.2016] - ggpubr [@Kassambara.2023] - openxlsx [@Schauberger.2023] - plyr [@Wickham.2011] - stringr [@Wickham.2023b] - tidyverse [@Wickham.2019].

Vor Beginn der Analysen müssen sowohl alle Packages (siehe @sec-Packages) geladen als auch die Working directory (siehe @sec-Working-Directory) gesetzt werden.

Der Code ist 'lauffähig', d.h. kann fehlerfrei gerechnet werden, so die vorgegebene Datenstruktur beibehalten wird.

Die verwendeten Daten wurden mit dem Bruker Tracer 5g (900F3980) des Departments für Kulturwissenschaften und Altertumskunde der Ludwig-Maximilians-Universität München im MudrockDual-Modus mit 300 Sekunden Messzeit (jeweils 150 Sekunden für Mudrock Major und Trace) und einem 8 mm Kollimator aufgenommen. Die Datenaufnahme erfolgte durch M. Schauer am 26. Januar 2023 (20°C; relative Luftfeuchte 40%). Da es sich um eine Pilotstudie handelt, erfolgte pro Messposition (Oberfläche, Auftrag, alter Bruch etc.) nur eine Messung. Zudem wurde keine Probenaufbereitung durchgeführt, weshalb die erzeugten Daten ausschließlich im qualitativen Sinne zum Ziehen von Vergleichen innerhalb der Messserie verwendet werden können: Sie bieten einen Einblick in das Material und können als Ausgangspunkt für weitere Analysen dienen.

Das hier abgebildete Verfahren der Datenkontrolle folgt grundsätzlich jenem nach Schauer [-@Schauer.2102023] bzw. Schauer/Amicone [-@Schauer.2024].Da im Rahmen der Untersuchung Einzelmessungen durchgeführt wurden, wird in diesem Falle jedoch die prozentuale relative Messunsicherheit (RU) als Prüfvariable berechnet, wobei der jeweilige Messwert sowie die 2δ-Messunsicherheit als Grundlage dienten. Die angewandte Formel lautet somit RU = (Messwert/Messunsicherheit)\*100 [vgl. auch @Harris.2007]. Dabei darf RU einen Wert von 20 % nicht überschreiten. Zudem darf die Überschreitung des Grenzwertes pro Probe und chemischem Element in maximal 20 % der Fälle erfolgen.

## Packages {#sec-Packages}

::: {.callout-note title="Hinweise" collapse="true" appearance="minimal"}
Nach <https://www.projectpro.io/recipes/load-package-r> and <https://www.dataquest.io/blog/install-package-r/>
:::

```{r}
# Laden der notwendigen R-packages
library(basictabler)
library(cowplot)
library(data.table)
library(stringr)
library(dplyr)
library(forcats)
library(ggplot2)
library(ggpubr)
library(openxlsx)
library(plyr)
library(tidyverse)
```

## Set working directory {#sec-Working-Directory}

::: {.callout-note title="Hinweise" collapse="true" appearance="minimal"}
Nach <https://martinctc.github.io/blog/rstudio-projects-and-working-directories-a-beginner%27s-guide/>, <https://www.grainge.org/pages/authoring/relative_paths/relative_paths.htm> und <https://r4ds.had.co.nz/workflow-projects.html>
:::

```{r}
knitr::opts_knit$set(root.dir = "./")
```

## Datenkontrolle

### Standardabweichung

#### Standardabweichung der Einzelmessungen

Messungen KF28_1, KF266_1 und KF266_2 sind auffällig häufig von zu hohen Messfehlern betroffen.

```{r}
# Laden der Daten
data<- read.csv("../Daten//Farchant_Bruker_Gesamt.csv")

# Definieren der zur Berechnung zu verwendenden Spalten
variables <- c("Si", "Ti", "Al", "Fe", "Mn", "Mg", "Ca", "K", "P", "S", "V", "Cr", "Co", "Ni", "Cu", "Zn", "As", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Ba", "Pb", "Th", "U")

# Identifizieren der Einzelmessungen deren Messwert niedriger ist als +- zugehöriger Messfehler (.Error)
for (variable in variables) {
  filter <- data[data[, variable] < 2*data[, paste0(variable, ".Error")], ]
  count <- nrow(filter)
  assign(paste0("data", variable), filter)
  assign(paste0("count", variable), count)
  print(paste0("count", variable, ": ", count))}

# Erstellen einer Liste relevanter Elemente (max. Anzahl auffällige Messungen = 4)
l = list(dataV,dataCr,dataCo,dataCu,dataNb,dataMo,dataPb,dataTh)
Tab1<-rbindlist(l, use.names=TRUE, fill=TRUE)
```

### Prozentuale relative Messunsicherheit (RU)

#### Berechnen der RU pro Probe

```{r}
# Laden der Daten
data1<- read.csv("../Daten//Farchant_Bruker_Gesamt.csv")

# Berechnen des Mittelwertes pro Probe (relevant für die dreifache Standardmessung)
Mean <- data1 %>%
  group_by(SAMPLE) %>%
  dplyr::summarize(across(everything(), ~mean(.x, na.rm = TRUE), .names = "{.col}"))

# Auswahl der zur Berechnung zu verwendenden Spalten
variables <- colnames(Mean)[17:71]

# Erstellen der Formel zur Berechnung der RU
calculate_RU <- function(variable) {
  mean_col <- variable
  sd_col <- paste0(variable, ".Error")
  if(mean_col %in% colnames(Mean) & sd_col %in% colnames(Mean)) {
    return((Mean[[sd_col]] / Mean[[mean_col]]) * 100)
  } else {
    return(rep(NA, nrow(Mean)))
  }
}

# Berechnen der RU für die zuvor definierten Variablen und anhängen dieser Spalten an die Datei Mean
for (variable in variables) {
  Mean[[paste0(variable, "_RU")]] <- calculate_RU(variable)
}

# Ersetzen der Informationen NA, INF und RU > 99 mit dem Wert 99
data_RU <- lapply(Mean, function(x) {
  x[is.na(x) |is.infinite(x) | x > 99] <- 99
  return(x)
})

# Definieren von data_RU als Tabelle
data_RU <- as.data.frame(data_RU)

# Definieren der Spalte SAMPLE der Datei Mean als Variable SAMPLE
SAMPLE<-Mean$SAMPLE

# Zusammenfügen von SAMPLE und data_RU in eine gemeinsame Tabelle
data2<-cbind(SAMPLE,data_RU)

# Auswahl der weiter zu analysierenden Elemente
data_RU<-data2[,c("SAMPLE","Si_RU","Ti_RU","Al_RU","Fe_RU","Mn_RU","Mg_RU","Ca_RU","K_RU","P_RU","S_RU","As_RU","V_RU","Cr_RU","Ni_RU","Cu_RU","Zn_RU","Rb_RU","Sr_RU","Y_RU","Zr_RU","Nb_RU","Ba_RU","Pb_RU")]
```

#### Erstellen zahlcodierte Tabelle

```{r}
# Ersetzen aller Werte <=20 mit '0' und >20 mit '1'
Prüfung_RU_0 <- cbind(apply(data_RU[,2:24], 2, function(x) ifelse(x <= 20, 0, x)), data_RU[,1])
Prüfung_RU_0_1 <- cbind(apply(Prüfung_RU_0[,1:23], 2, function(x) ifelse(x > 20, 1, x)), Prüfung_RU_0[,24])

# Definieren von Prüfung_RU_2 als Tabelle, Bennen der Spalte 24 als SAMPLE und definieren aller RU-Werte als Zahl
Prüfung_RU_2 <- as.data.frame(Prüfung_RU_0_1)
names(Prüfung_RU_2)[24] <- "SAMPLE"
Prüfung_RU_2[,1:23] <- sapply(Prüfung_RU_2[,1:23],as.numeric)

# Berechnen der Summe der zu hohen RU pro Zeile (Proben) unter Ausschluß von Elemente mit zu vielen zu hohen RUs
Prüfung_RU_3<-Prüfung_RU_2%>%mutate(Zeilensumme=rowSums(.[, !(names(.) %in% c("SAMPLE","Mg_RU","As_RU","V_RU","Cr_RU","Nb_RU","Ba_RU","Pb_RU"))]))

# Berechnen der Summe zu hohen RU pro Spalte (Elemente)
spaltensummen <- colSums(Prüfung_RU_3[, !colnames(Prüfung_RU_3) %in% "SAMPLE"])
spaltensummen <- c(spaltensummen,0)
Prüfung_RU_4 <- rbind(Prüfung_RU_3,spaltensummen)

# Berechnen des prozentualen Anteils der zu hohen RUs pro Spalte (Elemente) 
Zeilenanzahl<-(nrow(Prüfung_RU_4)-1)
Prozent<-(100/Zeilenanzahl*spaltensummen)
Prüfung_RU_5 <- rbind(Prüfung_RU_4,Prozent)

# Berechnen des prozentualen Anteils der zu hohen RUs pro Zeile (Proben) 
Zeilensumme<-Prüfung_RU_5$Zeilensumme
Spaltenanzahl<-ncol(Prüfung_RU_5)
Zeilenprozent<-(100/(Spaltenanzahl-9)*Zeilensumme)
Prüfung_RU_6 <- cbind(Prüfung_RU_5,Zeilenprozent)

# Runden aller Zahlenwerte auf ganze Zahlen und sortieren der Tabelle
ganze_Zahlen <- function(x) {if(is.numeric(x)) {return(round(x))} else {
    return(x)}}
Prüfung_RU_6 <- as.data.frame(lapply(Prüfung_RU_6,ganze_Zahlen))
Prüfung_RU_6<-Prüfung_RU_6[order(Prüfung_RU_6$Zeilenprozent, decreasing = TRUE),]

# Einfügen der Bezeichnung Spaltensumme und Spaltenprozent in die jeweilig zugehörige Zeile in Spalte SAMPLE
text<-"Spaltensumme"
Prüfung_RU_6[14,24]<-text

text<-"Spaltenprozent"
Prüfung_RU_6[15,24]<-text

# Ausgabe der Berechnungen
Prüfung_RU_6

# Export der Daten
write.csv(Prüfung_RU_6,"../Daten//Farchant_Bruker_RU_Zahlkodiert.csv")
```

#### Erstellen farbcodierte Tabelle

::: {.callout-note title="Hinweise" collapse="true" appearance="minimal"}
Nach <https://cran.r-project.org/web/packages/basictabler/vignettes/v04-styling.html> und <https://cran.r-project.org/web/packages/basictabler/vignettes/v07-excelexport.html>
:::

```{r}
# Erstellen einer Tabelle
tbl <- BasicTable$new()

# Befüllen mit den Daten aus data_RU
tbl$addData(data_RU, firstColumnAsRowHeaders=TRUE)

# Definieren der Zellen für welche eine Farbkodierung erfolgen soll
cells <- tbl$getCells(rowNumbers=2:14, columnNumbers=2:24, matchMode="combinations")

# Definieren und Anwenden des Farbkodes: rot wenn RU > 20, schwarz wenn RU <= 20
tbl$mapStyling(cells=cells, styleProperty="background-color", valueType="color",mapType="logic",mappings=list("v>20", "red", "v<=20", "black"))
tbl$renderTable()

# Erstellen einer exportfähigen .xlsx-Datei
wb <- createWorkbook(creator = Sys.getenv("M. Schauer"))
addWorksheet(wb, "Data")
tbl$writeToExcelWorksheet(wb=wb, wsName="Data",topRowNumber=1, leftMostColumnNumber=1, applyStyles=TRUE)

# Export der Datei
saveWorkbook(wb, file = "../Daten//Farchant_Bruker_RU_Farbkodiert.xlsx", overwrite = TRUE)
```

## Datenaufbereitung

Die MudrockDual-Spektren des Bruker Tracer 5g wurden zunächst mit der Software CalToolkitEC (v. 1.57.1) in die Spektren der einzelnen Modi MudrockMajor und -Trace aufgeteilt und dann in Artax (v. 8.0.0.476) weiterverarbeitet. Es erfolgte eine Peak-Identifikation zur Modellierung der Spektren mithilfe Bayesscher Dekonvolution. Darauf basierend wurde eine Berechnung der Nettozählrate, d.h. der Anzahl von Photonen pro chemischem Element nach Abzug von Elementinterferenzen und Hintergrund, durchgeführt. Zur Verbesserung der Vergleichbarkeit wurde abschließend eine Normierung auf das Targetmaterial, im Falle des Tracers 5g Rhodium, durchgeführt. Nur dieser letzte Schritt der Datenaufbereitung erfolgte in R. Der Prozess folgt [@Drake.2014].

### Normierung Bruker Artax Mudrock Major Daten

```{r}
# Laden der Daten
data<- read.csv("../Daten//Farchant_BrukerArtax_Major.csv")

# Erstellen der Funktion, um neue Spaltennamen zu generieren
Spaltentitel <- function(original_cols, suffix_col) {
  paste0(original_cols, "_div_", suffix_col)
}

# Berechnung der neuen Spalten
data1 <- data
new_cols <- lapply(6:36, function(i) {
  data[[i]] / data[[24]]
})

# Namen der neuen Spalten festlegen
new_colnames <- paste0("temp_", names(data)[6:36])
names(new_cols) <- new_colnames

# Hinzufügen der neuen Spalten zum DataFrame
data1 <- cbind(data1, new_cols)

# Neue Spaltennamen generieren
new_colnames <- Spaltentitel(names(data)[6:36], names(data)[24])

# Entfernen von '.K12' und '.L1' aus den Spaltentiteln und ersetzen von '_div_' mit '_'
new_colnames <- str_replace_all(new_colnames, "\\.(K12|L1)", "")
new_colnames <- str_replace_all(new_colnames, "_div_", "_")

# Anwenden der neuen Spaltennamen
colnames(data1)[which(startsWith(colnames(data1), "temp_"))] <- new_colnames

# Auswahl der relevanten Spalten aus data1
data2 <- data1[c(1:5, 37, 41, 44:45, 49, 53, 58, 61)]

# Export der Daten
write.csv(data2,"../Daten//Farchant_BrukerArtax_Major_norm.csv")
```

### Normierung Bruker Artax Mudrock Trace Daten

```{r}
# Laden der Daten
data<- read.csv("../Daten//Farchant_BrukerArtax_Trace.csv")

# Erstellen der Funktion, um neue Spaltennamen zu generieren
Spaltentitel <- function(original_cols, suffix_col) {
  paste0(original_cols, "_div_", suffix_col)
}

# Berechnung der neuen Spalten
data1 <- data
new_cols <- lapply(6:49, function(i) {
  data[[i]] / data[[33]]
})

# Namen der neuen Spalten festlegen
new_colnames <- paste0("temp_", names(data)[6:49])
names(new_cols) <- new_colnames

# Hinzufügen der neuen Spalten zum DataFrame
data1 <- cbind(data1, new_cols)

# Neue Spaltennamen generieren
new_colnames <- Spaltentitel(names(data)[6:49], names(data)[33])

# Entfernen von '.K12' und '.L1' aus den Spaltentiteln und ersetzen von '_div_' mit '_'
new_colnames <- str_replace_all(new_colnames, "\\.(K12|L1)", "")
new_colnames <- str_replace_all(new_colnames, "_div_", "_")

# Anwenden der neuen Spaltennamen
colnames(data1)[which(startsWith(colnames(data1), "temp_"))] <- new_colnames

# Auswahl der relevanten Spalten aus data1
data2 <- data1[c(1:5,62,73,75,85,89,91,92)]

# Export der Daten
write.csv(data2,"../Daten//Farchant_BrukerArtax_Trace_norm.csv")
```

## Datenauswertung - Histogramme

### Abb_Hauptelemente

```{r}
# Laden der Daten
data <- read.csv("../Daten//Farchant_BrukerArtax_Major_norm.csv")

# Definieren der einzigartigen Werte zur Sortierung
unique_values <- unique(data[[3]])

# Erstellen des Sortierungsvektors
Sortierung <- unique(c("alte Ofl. Helle Ware",
                        "alte Ofl. Roetliche Ware",
                        "alte Ofl. Wandputz",
                        "Ofl. Auftrag",
                        "Ofl. Grau",
                        "Ofl. Schwarz",
                        "Ofl. Weiss Helle Ware",
                        "Ofl. Weiss Wandputz"))

# Durchführen der Sortierung und setzen von unbekannten Werten auf NA
data1 <- data %>%
  mutate(custom_sort = factor(data[[3]], levels = Sortierung)) %>%
  arrange(custom_sort) %>%
  select(-custom_sort)

# Auswahl der darzustellenden Spalten
element_names <- colnames(data1)[c(7:14)]

# Definieren der Farben
meine_farben <- c("orange","firebrick","red","darkred","grey55","black","wheat","pink")

# Erstellen der Histogramme
plots <- list()

for(element in element_names) {
  sorted_data <- data1 %>%
    arrange(Messstelle)
  
  hist_plot <- ggplot(data = sorted_data, aes(x = fct_reorder(SAMPLE, Messstelle), y = !!sym(element), fill = Messstelle)) +
    geom_bar(stat = "identity", position = "dodge", color = "black") +
    labs(x = "Probe",
         y = element,
         fill = "Messstelle") +
    scale_fill_manual(values = meine_farben) +
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
  
  print(hist_plot)
  
  plots[[element]] <- hist_plot
}

# Erstellen der Legende
legend <- cowplot::get_legend(plots[[1]])

# Entfernen der Legende von jedem Plot und ersetzen durch NULL
for (i in seq_along(plots)) {
  plots[[i]] <- plots[[i]] + theme(legend.position = "none")
}

# Anordnen von Plots und Legende
multiplot <- cowplot::plot_grid(plotlist = plots, legend, ncol = 2)

# Exportieren der Abbildung
ggsave("../Abbildungen//Abb_Hauptelemente.eps", plot = multiplot, device = "eps", width = 15.3, height = 25, units = "cm")
```

### Abb_Spurenelemente

```{r}
# Laden der Daten
data <- read.csv("../Daten//Farchant_BrukerArtax_Trace_norm.csv")

# Definieren der einzigartigen Werte zur Sortierung
unique_values <- unique(data[[3]])

# Erstellen des Sortierungsvektors
Sortierung <- unique(c("alte Ofl. Helle Ware",
                        "alte Ofl. Roetliche Ware",
                        "alte Ofl. Wandputz",
                        "Ofl. Auftrag",
                        "Ofl. Grau",
                        "Ofl. Schwarz",
                        "Ofl. Weiss Helle Ware",
                        "Ofl. Weiss Wandputz"))

# Durchführen der Sortierung und setzen von unbekannten Werten auf NA
data1 <- data %>%
  mutate(custom_sort = factor(data[[3]], levels = Sortierung)) %>%
  arrange(custom_sort) %>%
  select(-custom_sort)

# Auswahl der darzustellenden Spalten
element_names <- colnames(data1)[c(7:13)]

# Definieren der Farben
meine_farben <- c("orange","firebrick","red","darkred","grey55","black","wheat","pink")

# Erstellen der Histogramme
plots <- list()

for(element in element_names) {
  sorted_data <- data1 %>%
    arrange(Messstelle)
  
  hist_plot <- ggplot(data = sorted_data, aes(x = fct_reorder(SAMPLE, Messstelle), y = !!sym(element), fill = Messstelle)) +
    geom_bar(stat = "identity", position = "dodge", color = "black") +
    labs(x = "Probe",
         y = element,
         fill = "Messstelle") +
    scale_fill_manual(values = meine_farben) +
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
  
  print(hist_plot)
  
  plots[[element]] <- hist_plot
}

# Erstellen der Legende
legend <- cowplot::get_legend(plots[[1]])

# Entfernen der Legende von jedem Plot und ersetzen durch NULL
for (i in seq_along(plots)) {
  plots[[i]] <- plots[[i]] + theme(legend.position = "none")
}

# Anordnen von Plots und Legende
multiplot <- cowplot::plot_grid(plotlist = plots, legend, ncol = 2)

# Exportieren der Abbildung
ggsave("../Abbildungen//Abb_Spurenelemente.eps", plot = multiplot, device = "eps", width = 15.3, height = 25, units = "cm")
```

## References

::: {#refs}
:::
