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

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(pXRF)

Feinkalibration

Um die Daten an die Messgenauigkeit einer stationären XRF anzunähern, bzw. um die spezifischen systematischen Abweichungen unseres speziellen Gerätes auszugleichen, kann eine Feinkalibration durchgeführt werden. Hierzu benötigt man die Messungen der gleichen Materie mittels beider Geräte. Die Messergebnisse können dann in Beziehung gebracht werden, und Korrekturfaktoren können erstellt werden.

Hierzu benötigen wir zuerst einen Datensatz, der beide Messungen für eine grössere Zahl von Objekten enthält. Diese Daten können unterschiedlich strukturiert sein, wir haben in unserem Fall eine Tabelle, die die Elementmessungen enthält, und zusätzlich einen Identifizierer für die Messmethode und einen für die Probe enthält.

fine_calibration_data_mining.raw_data <- read.csv("../inst/extdata/fine_calibration_data_mining.raw_data.csv")

Nehmen wir den Datensatz mal in zwei einzelne Variablen auseinander

pxrf <- fine_calibration_data_mining.raw_data[fine_calibration_data_mining.raw_data$method=="pXRF",]
wdrfa <- fine_calibration_data_mining.raw_data[fine_calibration_data_mining.raw_data$method=="WD RFA",]

Wenn wir jetzt auf der y-achse die Messungen des stationären, und auf der x-Achse die des pXRF Gerätes plotten (erstmal für ein Element), dann können wir uns diese Beziehung ansehen:

plot(pxrf$Si,
     wdrfa$Si,
     xlab = "pXRF",
     ylab = "WD RFA")

Betrachten wir nur den Datenbereich, so gibt es hier eine deutliche Streuung. Wenn wir den Blickwinkel erweitern, und auch die Nullstelle miteinbeziehen, dann sieht die Beziehung schon linearer aus:

plot(pxrf$Si,
     wdrfa$Si,
     xlab = "pXRF",
     ylab = "WD RFA",
     xlim = c(0,400000),
     ylim = c(0,400000),
     asp = 1
     )

Hieraus können wir nun verschiedene Dinge bestimmen. Zum einen können wir "Übersetzungsfaktoren" bestimmen, die angeben, wie die Messung auf dem einen Gerät in eine Messung auf dem anderen zu überführen wäre. Das sind die Koeffizienten. Das andere ist die Güte, mit der eine solche Übertragung möglich ist. Dies wird allgemein mit dem Determinationskoeffizienten bestimmt, der auch als R² abgekürzt wird.

Beginnen wir mit den Koeffizienten. Hierzu bestimmen wir die Linie, die mit dem geringsten Abstand zu allen Punkten hindurch führt. Dies nennt man dann auch Lineares Model (lm).

lineares_modell <- lm(wdrfa$Si ~ pxrf$Si)
plot(pxrf$Si,
     wdrfa$Si,
     xlab = "pXRF",
     ylab = "WD RFA",
     xlim = c(0,400000),
     ylim = c(0,400000),
     asp = 1
     )
abline(lineares_modell)

Aus diesem Modell können wir nun die Übersetzungsparameter herleiten:

lineares_modell$coefficients

Dieses Modell besagt, dass wenn das pXRF 0 messen würde, dies mit dem Wert 10820 im stationären XRF zu übersetzen wäre. Und für jede weitere Einheit der Messung auf dem pXRF sind 0.67 auf diesen Wert aufzuaddieren.

Dies ergibt sich aber unter anderem dadurch, dass wir keine Messungen mit so geringen Silizium-Werten vorliegen haben. Wenn wir daher davon ausgehen, das beim Nichtvorhandensein des Elementes beide Geräte 0 anzeigen, dann ergibt sich ein etwas anderes Modell:

lineares_modell <- lm(wdrfa$Si ~ 0 + pxrf$Si)
plot(pxrf$Si,
     wdrfa$Si,
     xlab = "pXRF",
     ylab = "WD RFA",
     xlim = c(0,400000),
     ylim = c(0,400000),
     asp = 1
     )
abline(lineares_modell, col="red")

Hier geben wir an, dass wir davon ausgehen, dass es keinen Versatz für den 0-Wert gibt (daher 0 in der Formel), und dass alle Effekte nur über einen systematischen Fehler des pXRF zu erklären sind.

Daraus ergeben sich folgende(r) Koeffizient:

lineares_modell$coefficients

Das bedeutet also, dass wir den Wert des pXRF für Silizium mal 1.063579 nehmen, um den equivalenten Wert für das WD-XRF zu erhalten, der mit diesem gemessen worden wäre.

Für dieses Modell können wir (ebenso wie für das vorgehende Modell) den Determinationskoeffizienten berechnen (lassen):

summary(lineares_modell)$r.squared

Das bedeutet, dass sich nahezu 99% der Variabilität in den WD-XRF-Daten durch die pXRF Daten erklären lassen, und nur 1% zufällige Streuung ist.

Outlier

Schauen wir uns das ganze mal für Kupfer an:

plot(pxrf$Cu,
     wdrfa$Cu,
     xlab = "pXRF",
     ylab = "WD RFA"
     )

Hier gibt es zwei Punkte, die ganz klar und deutlich sich von allen anderen absetzen. Das können Daten- oder Messfehler sein. Diese Extremwerte beeinflussen sehr stark unser Lineares Modell, und machen es daher schlecht angepasst für den grössten Teil der Daten:

plot(pxrf$Cu,
     wdrfa$Cu,
     xlab = "pXRF",
     ylab = "WD RFA"
     )
abline(lm(wdrfa$Cu~pxrf$Cu))

Es kann also sehr sinnvoll sein, diese Extremwerte nicht in die Berechnung einzubeziehen. Man kann nun (und das ist sicher die bessere Methode) diese Werte einzeln bestimmen. Oder man überlässt dies einem Algorithmus.

Ein Ansatz ist der, der auch im Boxplot verwendet wird:

boxplot(wdrfa$Cu)

Die meisten Messwerte sind im unteren Bereich, aber zwei Werte sind deutlich höher gemessen im WDXRF, was auf eine Verunreinigung hinweisen könnte. Das Mass ist hier die Box des Boxplots: Diese entspricht den werten der mittleren Hälfte der Daten. Wenn ein Wert mehr als 1.5x der Grösse dieser Spanne, gemessen von derer oberen bzw. unteren Grenze, abweicht, so wird er als Ausreisser verstanden.

Wir könnten die in dem Boxplot eingebaute Funktion nutzen, oder wir schreiben sie zum besseren Verständnis kurz selbst:

iqr <- IQR(wdrfa$Cu) # interquartile
quantiles <- quantile(wdrfa$Cu)
outlier <- wdrfa$Cu > quantiles[4] + iqr * 1.5 | wdrfa$Cu < quantiles[2] - iqr * 1.5
wdrfa$Cu[outlier]

Das identifiziert uns die beiden Outlier. Jetzt verbinden wir diesen Ansatz in einer Funktion mit dem Linearen Modell zu einer robusten Regression:

outlier <- function(x) {
iqr <- IQR(x) # interquartile
quantiles <- quantile(x)
outlier <- x > quantiles[4] + iqr * 1.5 | x < quantiles[2] - iqr * 1.5
outlier
}

robust_lm <- function(x,y) {
  outlier_x <- outlier(x)
  outlier_y <- outlier(y)
  lm(y[!(outlier_x | outlier_y)] ~ 0 + x[!(outlier_x | outlier_y)])
}

robust_lm(pxrf$Cu, wdrfa$Cu)

Und nun noch darstellen, erst gegen den gesamten Datensatz:

plot(pxrf$Cu,
     wdrfa$Cu,
     xlab = "pXRF",
     ylab = "WD RFA"
     )
abline(robust_lm(pxrf$Cu, wdrfa$Cu))

Und nun gegen den Datensatz, der keine Outlier mehr enthält:

plot_without_outlier <- function(x,y,...) {
  outlier_x <- outlier(x)
  outlier_y <- outlier(y)
  plot(x[!(outlier_x | outlier_y)],y[!(outlier_x | outlier_y)],...)
}

plot_without_outlier(pxrf$Cu,
     wdrfa$Cu,
     xlab = "pXRF",
     ylab = "WD RFA"
     )
abline(robust_lm(pxrf$Cu, wdrfa$Cu))

residuals(lineares_modell)

Final

Mit ein bisschen r-Magie können wir das nun effizient für alle Elemente berechnen lassen:

correction <- sapply(3:ncol(pxrf), function(x)
  robust_lm(pxrf[,x],wdrfa[,x])$coefficients
  )
names(correction) <- colnames(pxrf)[3:ncol(pxrf)]
correction

Und gleich noch die Güte der Übertragungen als R²

r_squared <- sapply(3:ncol(pxrf), function(x)
  summary(robust_lm(pxrf[,x],wdrfa[,x]))$r.squared
  )
names(r_squared) <- colnames(pxrf)[3:ncol(pxrf)]
sort(r_squared)