In MartinHinz/pXRF: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_knit$set(root.dir = '.')

library(pXRF)

Zuerst laden wir die bereinigten Daten ein:

# Durch interaktive Wahl des Benutzers
#source_csv <- file.choose()

# Oder durch einen fest angestellten Pfad
source_csv <- "data_feinkalibriert.csv"

all_data <- read.csv(source_csv, row.names = "SAMPLE")

An dieser Stelle gehe ich davon aus, das eindeutige Probennamen existieren. Dies sollte auf jeden Fall so sein, und sollte das Kommando zu Laden einen Fehler auswerfen, dann müsste hier geschaut werde!

Gegebenfalls können wir die Daten nun noch auf ein Kriterium hin filtern, z.B. bezüglich des Messmodus (Type). Die zu filternde Spalte und den zu filternden Wert kann man entsprechend den Notwendigkeiten austauschen.

data_ausw<-subset(all_data, Type == "Mining")

Als nächstes wollen wir die Daten hinsichtlich der Elementzusammensetzungen in Streudiagrammen darstellen, wobei immer die Werte bezüglich zweier Elemente gegeneinander geplottet werden.

Dazu kann man auf bestimmte Spalten, und damit auf bestimmte Elemente auswählen:

elements_selected <- c("Nb", "Zr", "Y", "Sr", "Rb",  "Zn",  "Fe",  "Cr", "V", "Ti",  "Ca", "K", "Al",  "Si", "Th", "Pb", "Cu", "Ni", "Mn", "Ba", "P","S", "Mg")

# Alle Elemente: "Nb", "Zr", "Y", "Sr", "Rb", "Th", "Pb", "Zn", "Cu", "Ni", "Fe", "Mn", "Cr", "V", "Ti", "Ba", "Ca", "K", "Al", "P", "Si", "S", "Mg"
# Th, Pb, Cu, Ni, Mn, Ba, und Mg oft unkonstant gemessen resp. P und S Kontamination, daher eher weglassen...

data <- data_ausw[,elements_selected]

Jetzt kann man direkt den Plot vornehmen:

pairs(data, main="Streudiagramme der Messungen")

Hier werden jetzt sehr viele Elemente gegeneinander geplotted. Man kann die Auswahl auch beschränken, durch Angabe einzelner Elemente, oder durch die Spaltennummern

pairs(data[,c("Ca", "K", "Al")], main="Streudiagramme der Messungen")

Falls man untergruppen deutlich machen möchte, kann man dies mittels Symbolen und/oder Farben tun. Ich habe die Spalte "NOTE" in der ursprünglichen Tabelle genutzt, um eine Gruppenzugehörigkeit mittels einzelner Werte zuzuordnen. Bei mir sind dies die einzelnen Materialgruppen. Dies kann aber auch alles andere mögliche sein, dass eine Gruppenzugehörigkeit für den jeweiligen Datensatz sinnvoll erscheinen lässt.

data_ausw$NOTE

Diese Informationen übertrage ich in eine eigene Variable, die ich groups nenne. Farbdarstellung bedingt dies als numerische Werte, daher lege ich eine weitere Variable groups_num an, in der den Gruppen einfach aufsteigende Zahlen zugeordnet werden. Schliesslich brauch ich für die automatische Farbdarstellung noch die gesamtzahl von Gruppen, die ich in groups_count notieren lasse.

groups <- factor(data_ausw$NOTE)
groups_num <- as.numeric(groups)
groups_count <- length(unique(groups))

Im folgenden benutze ich die Funktion rainbow, um mir automatisch eine bestimmte Anzahl (gleich der Anzahl der Gruppen) von Farben ausgeben zu lassen, die ich für das Einfärben verwenden kann. Gleichzeitig benutze ich die option "pch", um unterschiedliche Symbole für die einzelnen Gruppen darzustellen. Eigentlich braucht man nur eine der beiden Optionen, um eine Zuordnung im Plot und für die Auswertung treffen zu können...

pairs(data[,c("Ca", "K", "Al")],
      main="Streudiagramme der Messungen",
      col = rainbow(groups_count)[groups_num],
      pch = groups_num
      )

Weis man um die Zahl der Gruppen, und möchte man diese immer gleich färben, so bietet es sich an, die Farbpalette einmal zu definieren, um sie immer wieder verwenden zu können. Dies ist ganz einfach, die meisten gängigen Farbennamen erkennt R.

my_colors <- c("blue", "red", "green", "black")

pairs(data[,c("Ca", "K", "Al")],
      main="Streudiagramme der Messungen",
      col = my_colors[groups_num],
      pch = groups_num
      )

Alle bekannten Farbnamen kann man mit folgendem Befehl einsehen (den ich allerdings auskommentiert habe):

# colors()

Um die Zuordnung von Farben oder Symbolen im Plot anzuzeigen, kann man noch eine Legende hinzufügen:

pairs(data[,c("Ca", "K", "Al")],
      main="Streudiagramme der Messungen",
      col = rainbow(groups_count)[groups_num],
      pch = groups_num
      )
legend("topright", legend=unique(groups), pch = unique(groups_num), col = rainbow(groups_count)[unique(groups_num)])

Die Legende sitzt nicht schön (bzw. ausserhalb des Plotbereiches), daher platzieren wir sie unterhalb des Plots

pairs(data[,c("Ca", "K", "Al")],
      main="Streudiagramme der Messungen",
      col = rainbow(groups_count)[groups_num],
      pch = groups_num, oma=c(10,3,3,3)
      )

legend("bottom", legend=unique(groups), pch = unique(groups_num), col = rainbow(groups_count)[unique(groups_num)],
       xpd=T, horiz = T)

Das gleiche nochmal etwas schöner, unter Nutzung einer vordefinierten Funktion

library(car)
scatterplotMatrix(data[,c("Ca", "K", "Al")],
                  diagonal=list(method ="histogram", breaks="FD"),
                  smooth=FALSE)

Wir können auch zwei Elemente separat gegeneinander plotten:

plot(data$Ca,
     data$K)

Oder schöner:

scatterplot(data$Ca,
     data$K,groups=groups)

Mittelwerte und Standardabweichungen

Die Frage, ob sich die einzelnen Gruppen unterscheiden, macht sich ja daran fest, ob sich die Messwerte zwischen diesen unterscheiden. Dabei streuen die Messwerte innerhalb der Gruppen je Einzelmessung (gleiches oder unterschiedliches Objekt) ja auch mess- und materialbedingt. Das Kriterium, ob man Gruppen identifizieren kann, ist dabei, ob die Unterschiede innerhalb der Gruppen kleiner sind als die Unterschiede zwischen den Gruppen.

Hierzu sind zwei deskriptive Eigenschaften von Bedeutung: - Der Mittelwert (oder einer der Varianten von Mittelwerten) gibt an, welche Werte die Gruppe insgesamt aufweist - Die Standardabweichung gibt an, wie die Daten um diesen Mittelwert streuen.

Wir können uns diese Werte recht einfach mittels einer bequemen R-Funktion anzeigen lassen

aggregate(data, list(groups), mean)
aggregate(data, list(groups), median)
aggregate(data, list(groups), sd)

Dies bleibt jedoch in der puren Zahlenansicht nicht sehr intuitiv. Eine Graphische Darstellung wie ein Boxplot ist hier hilfreicher. Für einzelne Element geht das gut mit der Standard-Graphik

boxplot(data$Si ~ groups)

Für komplexere Darstellungen (viele Elemente + Gruppen gleichzeitig) gibt es die ggplot Bibliothek. Als erstes müssen wir dazu unsere Daten in das "Lange" format übertragen

library(reshape2)
data_long <- melt(cbind(data,groups))

Jetzt können wir diese dann auch mit Standard-Graphik darstellen, das bleibt aber unbefriedigend, da es sehr viele Elemente hat, die sich auch in ihren Messbereichen deutlich unterscheiden:

boxplot(value ~ variable + groups, data = data_long)

Besser ist es, wenn wir die einzelnen Plots per Element unterteilen

library(ggplot2)

ggplot(data = data_long) + geom_boxplot(aes(fill=groups, y = value)) + facet_wrap(.~variable, scales = "free_y")

Wollen wir nun schauen, ob sich 2 Gruppen in Bezug auf ihre Ausprägungen in Einzelnen Elementen signifikant unterscheiden, so können wir einen einfachen nichtparametrischen Test anwenden, wie z.B. den Wilcoxon Rang-Summen-Test.

wilcox.test(data$Si[groups=="Keramik"], data$Si[groups=="Silex"])

oder für alle Gruppen zugleich

pairwise.wilcox.test(data$Si, groups)

PCA

Um latente Variablen zu identifizieren, die mehrere (oder alle) Elemente gleichzeitig beeinflussen, bietet sich die Hauptkomponentenanalyse (principal component analysis, pca) an. Hierzu dürfen keine NA-Werte in den Daten vorhanden sein. Filtern wir sicherheitshalber erst mal dazu in dieser Hinsicht:

# 1 for rows, 2 for columns
na_columns <- apply(data,2,anyNA)

data.for_pca <- data[,!na_columns]

Jetzt können wir die Hauptkomponentenanalyse selbst durchführen:

data.pca <- prcomp(data.for_pca)

und darstellen

biplot(data.pca)

Häufig unterscheiden sich die einzelnen Elemente hinsichtlich der Grössenordnung ihrer Messwerte. Um besonders häufig enthaltene Element nicht überzubewerten, bzw. um Spurenelement nicht unterzubewerten, bietet es sich an, diese auf ein gleiches Mass zu normieren. Hierzu wird die sogenannte z-Transformation angewendet, die für alle Werte:

• den Mittelwert pro Element bestimmt und diesen von allen Elementen subtrahiert. Dadurch wird der Mittelwert für alle auf Null gesetzt.
• die Standardabweichung für alle Elemente bestimmt und durch diese teilt. Damit haben alle Elemente eine Standardabweichung von 1.

Dadurch werden die Grössenordnungen aller Elemente angeglichen, und ihr Einfluss auf die PCA damit ebenso. In R können wir das durch eine einfache Zugabe einer weiteren Option erreichen:

data.pca_scaled <- prcomp(data.for_pca, scale. = T)

Und wiederum dargestellt:

biplot(data.pca_scaled)

Eine weitere Möglichkeit, die Daten vorzubehandeln, um den Einfluss extremer Werte und Wertunterschiede zu verringern, ist die die Transformation, z.B. mittels des Logarithmus (zur Basis 10). Hierbei werden die ordinalen Unterschiede zwischen den einzelnen Messwerten bzw. Elementen nicht gänzlich aufgelöst, sondern nur abgeschwächt. Analog zur z-Transformation ist die Durchführung nicht kompliziert. Allerdings ist zu beachten, dass sich kein log10 von 0 bilden lässt. Daher bietet es sich an, auf jeden Werte einen (sehr kleinen) Wert aufzuaddieren, damit diese Problem umgangen wird.

data.pca_log10 <- prcomp(log10(data.for_pca+0.1))
biplot(data.pca_log10)

Schönere Darstellung mit ggplot

Es gibt verschiedene Pakete, mit denen die Durchführung einer PCA erleichtert, und/oder mit der die Darstellung schöner gestaltet werden kann. Ein Quasi-Standart im Moment ist die Graphikbibliothek ggplot2 und daraus abgeleitete Darstellungsformen. Diese Funktionalitäten kann man mittels Zusatz-Paketen installieren.

Zusatz-Pakete sind entweder in den offiziellen Paketquellen verfügbar, oder sie haben es noch nicht in diese geschafft, und müssen z.B. von GitHub installiert werden.

Ein Paket, dass speziell für die Darstellung von Hauptkomponenten, oder ähnlichen Verfahren, enwickelt ist, ist ggord. Dieses ist nicht aus den Offiziellen Quellen installierbar, sondern muss mittels eines weiteren Pakets, "devtools", von GitHub installiert werden.

Die folgenden Befehle für die Paket-Installation sind auskommentiert, damit sie nicht im Laufe des Scriptes hier jedes mal ausgeführt werden. Wenn gewünscht, bitte den Bereich hinter dem '#' markieren und auf Run klicken (oder anderweitig ausführen).

Als erstes installieren wir 'devtools':

# install.packages("devtools")

Jetzt machen wir das frisch installierte Paket devtools verfügbar, mittels des Befehls "library()":

# library(devtools)

Nun können wir das Paket ggord aus dem Repository auf GitHub (https://github.com/fawda123/ggord/) installieren:

# install_github('fawda123/ggord')

Alles weitere wieder unkommentiert, denn nun gehe ich davon aus, das ggord installiert ist. Wir laden ggord und führen die Darstellung durch. Zu bemerken sind folgende Parameter:

• der erste Wert in der Klammer gibt das darzustellende Resultat der Hauptkomponentenanalyse an
• der zweite Wert (vec_ext) gibt an, dass die Pfeile grösser (oder kleiner) dargestellt werden. 1 ist unverändert, werte über 1 vergrössern, werte unter 1 verkleinern
• der dritte Wert (exp) bewirkt, dass das Diagramm etwas grösser dargestellt wird, als der eigentliche Wertebereich ist. Das verhindert, dass Punkte zu dicht am Rand sind.

library(ggord)

ggord(data.pca_scaled, vec_ext = 5, exp = c(.1,.1))

Wir können auch die Probenamen verwenden:

ggord(data.pca_scaled, vec_ext = 5, exp = c(.1,.1), obslab = T)

Wenn wir jetzt vordefinierte Gruppen haben, können wir diese im Plot kenntlich machen

ggord(data.pca_scaled, grp_in = groups, vec_ext = 5, exp = c(.1,.1))

Ebenso können wir alle skalierten und transformierten PCAs darstellen:

ggord(data.pca_log10, groups, vec_ext = 5, exp = c(.1,.1))

Clusteranalyse

Die Clusteranalyse ist ein Werkzeug, um Gruppenzugehörigkeit aufgrund von unterschiedlichen Merkmalen mittels eines Algorithmus zu definieren. Hierbei werden aufgrund von unterschiedlichen Verfahren Ähnlichkeiten zwischen den einzelnen Objekten festgelegt, und dann die jenigen Objekte zu Gruppen zusammen gefügt, die sich am ähnlichsten sind.

Hierzu gibt es verschiedene Herangehensweisen. Eine Variante ist die so genannte hierarchische Clusteranalyse, bei der schrittweise Objekte in absteigender Ähnlichkeit zu einander zusammen gepackt werden. Hieraus entsteht ein so genannter Baum, der angebt, in welcher Reihenfolge Gruppierungen vorgenommen worden sind.

Da wir es mit Messwerten zu tun haben, können wir für die Bestimmung von Ähnlichkeiten die euklidische Distanz benutzen. Das ist das Distanzmass, dass man auch aus dem Alltag kennt. Man kann es als Abstand in Bezug auf X und Y Koordinaten sehen. Mathematisch können aber auch beliebig viele Dimensionen berücksichtigt werden.

Die Funktion, um diese Distanzen zu berechnen, ist die Funktion dist(). Ohne weitere Angabe der Distanzfunktion wird euklidische Distanz hiermit berechnet.

data.dist <- dist(data)

Die eigentliche Gruppierung in Cluster erfolgt nun mittels der Funktion hclust(). Auch für diese gibt es verschiedene Methoden. Für kritische Distanzen bietet sich die so genannte ward-Methode an.

data.hclust <- hclust(data.dist, method = "ward.D2")

Zur Darstellung des Clusterbaums kann man nun das entstandene Objekt einfach mit der Plotfunktion darstellen.

plot(data.hclust)

An den einzelnen Ästen des Baumes sind am Ende jeweils die Zeilen Namen der Objekte angetragen. Man kann diese jedoch auch zum Beispiel durch die Gruppenzugehörigkeit ersetzen. Hierzu wird die gewünschte Variable als "labels" mit angegeben

plot(data.hclust,labels = groups)

Eine hierarchische Clusteranalyse produziert nicht ausschliesslich eine Lösung, sondern einen Ähnlichkeitsbaum. Möchten wir uns nun auf eine Anzahl von Clustern festlegen, so müssen wir "diesen Baum fällen". Der Befehl hier zu lautet cutree().

data.clusters <- cutree(data.hclust, k=4)
data.clusters

Die daraus abgeleiteten Gruppen Zugehörigkeiten können wir nun nutzen, um eine Darstellung, zum Beispiel eine Hauptkomponentenanalyse, mit den abgeleiteten Gruppen ein zu färben.

plot(data.pca$x[,1],data.pca$x[,2], col = data.clusters)

Ein weiteres Clusterverfahren ist das so genannte kmeans-Verfahren. Dieses produziert häufig bessere Ergebnisse, verlangt aber, dass die Daten in euklidischen Distanz vorliegen. Zudem muss man vorher angeben, wie viel Cluster man erwartet. Genau diese Anzahl von Clustern wird dann auch eingeteilt.

data.kmeans <- kmeans(data, 4)

Auch diese Lösung können wir dann verwenden, um eine Darstellung entsprechend der Clusterzugehörigkeit ein zu färben.

plot(data.pca$x[,1],data.pca$x[,2], col = data.kmeans$cluster)

Diskriminanzanalyse (lda)

Die Diskriminanzanalyse ist ein Verfahren, um herauszufinden, welche Elemente zwei Gruppen in einem Datensatz am besten unterscheiden. Das heisst, wir haben schon eine Gruppeneinteilung, und möchten nun herausfinden, worin die Unterschiede zwischen diesen Gruppen bestehen.

library(MASS)
data.for_lda <- data
data.for_lda$groups <- groups
data.lda <- lda(groups ~ ., data = data.for_lda)
data.lda

Die Werte zeigen zuerst die Anteile der jeweiligen Gruppen, dann deren Mittelwerte bezüglich der Elemente, schliesslich die Diskriminanzkoeffizienten. Diese zeigen auf, ob ein Element eher auf die Zugehörigkeit zur einen oder zur anderen Gruppe hinweist. Deutlicher wird dies, wenn wir und das dazu noch graphisch darstellen:

plot(data.lda)

Schöner mit Paketen, die auf ggplot basieren:

ggord(data.lda, data.for_lda$groups, vec_ext = 40, exp = c(.1,.1))

MartinHinz/pXRF documentation built on Dec. 17, 2021, 3:15 a.m.