Einleitung

In den letzten Jahren erfreut sich R, insbesondere bei Unternehmen, immer größerer Beliebtheit. Während viele Nutzer bereits heute komplexe Probleme schnell und effizient mit Hilfe von R lösen können, steckt die adequate Präsentation der Ergebnisse vieler Orts noch in den Kinderschuhen. Zumeist werden Ergebnisse nach der eigentliche Analyse einfach in Word Dokumente kopiert und die entsprechenden Texte und Tabellen dort editiert. Dieser Workflow birgt einige offensichtliche Nachteile. Ändern sich Daten oder Modelle auch nur geringfügig, müssen entsprechende Texte und Tabellen umständlich geändert werden. Besonders im Zeitalter von Big Data und kurzlebigen Daten ist das ein nicht zu verachtender Nachteil.

Bereits seit einiger Zeit bietet R-Studio mit der Integration von knitr hierfür eine entsprechende Lösung an. knitr ermöglicht automatisierte Berichtserstellung direkt in R-Studio. Dazu werden R-Code und Text direkt innerhalb einer Datei geschrieben und später evaluiert. knitr unterstützt den Export in HTML, PDF und Markdown. Die Vorteile liegen auf der Hand:

  1. Die Verknüpfung von Analyse und Report ermöglicht es dem Nutzer auch bei größeren Änderungen der Daten die Berichte automatisiert anpassen zu können. Alle Grafiken, Tabellen und verknüpften Ergebnisse werden automatisch aktualisiert.
  2. Die vielfältigen Formate für den Export erlauben ein breites Anwendungsgebiet. So können zum Beispiel Berichte direkt ins Web gestellt (HTML) oder ausgedruckt (PDF) werden.
  3. Im Gegensatz zu Word Dateien bieten sowohl HTML als auch PDF eine standartisierte Darstellung auf allen Endgeräten. Nicht zu verachten ist auch die Kompatibilität beider Formate mit mobilen Devices wie Smartphones und Tables.

Reporting Automatisieren mit R

Ziel einer effizienten Automatisierung muss sein, das der Nutzer selbst bei größeren Änderungen der verwendeten Daten oder Modelle nur noch minimal Anpassungen von Hand erledigen muss. Um das zu erreichen werden sämliche relevanten Zahlen dynamisch in R berechnet, in Objekte abgelegt und auf diese Objekte dann im Fließtext zugegriffen.

Folgendes Beispiel soll dies illustrieren:

# Wir speichern zunächst unsere Ausgngsvariablen ab
x <- 1
y <- 1
# Im nächsten Schritt berechnen wir dann unsere Zielgröße und speichern auch diese in einem eigenen Objekt
z <- x + y

Wird \(z\) nun im Fließtext referenziert, wird der oben berechnete Wert automatisch ausgegeben.

Whisky, Whisky, Whisky

Wir bei STATWORX lieben natürlich Statistik, wir genießen aber auch gerne mal einen guten schottischen Single Malt Whisky. Inspiriert durch den Beitrag der Kollegen von Revolution Analytics, verbinden wir nun diese beiden Leidenschaften und stellen euch hier anhand eines Beispiels das Reporting mit knitr vor.

Für unsere Auswertung benötigen wir eine Reihe von Packages, die wir als erstes laden.

library(ggplot2)
library(xtable)
library(reshape2)
library(ggmap)
library(rgdal)
library(maptools)
library(RCurl)
library(gridExtra)
library(plyr)

Im nächsten Schritt laden wir den benötigten Datensatz von diesen Homepage.

# Daten von Homepage laden
temporaryFile <- tempfile()
download.file("https://www.mathstat.strath.ac.uk/outreach/nessie/datasets/whiskies.txt", destfile = temporaryFile, method="curl")
data <- read.table(temporaryFile, sep = ",", header = TRUE)
# Eventuell fehlende Werte ausschließen
data <- na.omit(data)
# Anzahl der Whiskies speichern
n.whiskies <- nrow(data)

In unserem Datensatz befinden sich nun insgesamt 86 verschiedene Whiskies (streng genommen sind es 86 Destillerien, jedoch befindet sich pro Destillerie jeweils nur ein Whisky im Datensatz). Whisky Kenner mögen uns dahingehend gändig sein, dass auch eine gewisse geschmackliche Varianz innerhalb der Destillerie durchaus üblich ist. Die Anzahl der Whiskies haben wir zuvor in dem Objekt n.whiskies gespeichert und können wir nun an jeder beliegen Stelle des Textes durch folgenden Code abrufen:

# `r n.whiskies` (ohne #)

Werden nun beispielsweise weitere Whiskies zum Datensatz hinzugefügt, so ändert sich auch der Bericht automatisch sobald der R-Code erneut ausgeführt wird. In Verbindung mit dem automatischen Download der Daten zu Beginn, haben wir somit immer den aktuellsten Datensatz für unsere Analyse zur Verfügung. Ein kurzer Blick auf den Datensatz zeigt uns, dass für jeden Whisky neben Name und geografischer Position (Längen- und Breitengrade) auch eine Reihe von Geschmacksindikatoren aufgeführt sind. Diese sind auf einer Skala von 0 bis 4 codiert, wobei 4 die höchste Ausprägung der jeweiligen geschmacklichen Dimension darstellt.

# Die ersten Zeilen des Datensatzes ausgeben
head(data)
##   RowID  Distillery Body Sweetness Smoky Medicinal Tobacco Honey Spicy
## 1     1   Aberfeldy    2         2     2         0       0     2     1
## 2     2    Aberlour    3         3     1         0       0     4     3
## 3     3      AnCnoc    1         3     2         0       0     2     0
## 4     4      Ardbeg    4         1     4         4       0     0     2
## 5     5     Ardmore    2         2     2         0       0     1     1
## 6     6 ArranIsleOf    2         3     1         1       0     1     1
##   Winey Nutty Malty Fruity Floral    Postcode Latitude Longitude
## 1     2     2     2      2      2  \tPH15 2EB   286580    749680
## 2     2     2     3      3      2  \tAB38 9PJ   326340    842570
## 3     0     2     2      3      2   \tAB5 5LI   352960    839320
## 4     0     1     2      1      0  \tPA42 7EB   141560    646220
## 5     1     2     3      1      1  \tAB54 4NH   355350    829140
## 6     1     0     1      1      2    KA27 8HJ   194050    649950

Geschmacks-Clustering mit k-Means

Ziel der folgenden Analyse soll sein, die verschiedenen Whiskies anhand deren Geschmacksprofil in ähnliche Gruppen (Cluster) einzuteilen. Hierzu verwenden wir den k-Means Cluster-Algorithmus, der basierend auf der Distanz der einzelnen Whiskies zu den sog. Cluster-Zentroiden ähnliche Geschmacksprofile in homogene Gruppen zusammenfasst. Für interessierte Leser finden hier mehr Informationen zum k-means Verfahren. Für das Clustern nach k-means müssen zunächst zwei Schritte durchgeführt werden:

  1. Standardisieren der zum Clustern verwendeten Eingangsgrößen.
  2. Bestimmung der Anzahl der zu verwendenden Cluster.

Da die relevanten Variablen in unserem Datensatz bereits alle auf der selben Skala gemessen wurden, können wir direkt zu der Bestimmung der Cluster-Anzahl übergehen. Eine Möglichkeit, die optimale Anzahl der Cluster zu bestimmen ist sich die Sum of Squares (Streuung) innerhalb der Gruppen für eine beliegige Anzahl Cluster zu berechnen. In R kann dies durch die folgenden Befehle erreicht werden:

# Berechnung der Sum of Squares und abspeichern dieser in einem Data Frame
wss <- (n.whiskies - 1) * sum(apply(data[, 3:14], 2, var))
for (i in 2:15) {
  wss[i] <- sum(kmeans(data[, 3:14], centers = i)$withinss)
  }
# Abspeichern in Data Frame
df <- data.frame(wss = wss,
                 number = seq(1,15, 1),
                 diff = rep(NA, 15),
                 percent = rep(NA, 15)
                 )

Um nun die genaue Anzahl der Cluster zu bestimmen berechnen für jeden potenzieller Cluster-Lösung die prozentuelle Änderung der Innergruppen-Streuung zur vorherigen Cluster-Lösung. Dort wo durch einen zusätzlichen Cluster die Abnahme der Innergruppen-Streuung am geringsten ist (“Knick” in der Kurve), findet sich die optimale Anzahl der zu berechnenden Cluster.

# Bestimmung der Differenz zum vorherigen Cluster
for (i in 2:15) {
  df$diff[i] <- df$wss[i-1] - df$wss[i]
  df$percent[i] <- df$diff[i] / df$wss[1]
}
# Grenzwert definieren
threshold <- 0.05
vector.cluster <- df$number[df$percent > threshold]
n.cluster <- vector.cluster [length(vector.cluster)]

Wir definieren unser Schwelle nach der wir keine neuen Cluster mehr zulassen mit 5 Prozent und erhalten somit 4 Cluster als optimale Cluster-Anzahl. Durch diese Schritte wird die Bestimmung der Cluster automatisch nach den von uns vorgegeben Parametern ausgeführt und uns reportet. Diese Schritte lassen sich auch mit Hilfe von ggplot2, einem R-Package zur Datenvisualisierung, sehr schön darstellen.

# Wir erstellen einen Plot mit ggplot2 und fügen die vorher berechneten Sum of Squares, die Anzahl der Cluster sowie unseren "Threshold" ein
ggplot(df, aes(x = number, y = wss)) + geom_point() + geom_line() + xlab("Number of Clusters") + ylab("Within Groups um Sum of Squares") + geom_vline(x = n.cluster, color = "red", linetype = 4, alpha = .75) + theme_bw()

plot of chunk unnamed-chunk-8

Wie aus dem Plot hervorgeht, nimmt die Sum of Squares nach dem 4. Cluster deutlich weniger ab als zuvor. Somit entscheiden wir uns dafür, die Analyse mit 4 Whisky-Clustern durchzuführen.

# K-means Clustering mit 4 Whisky-Clustern
fit.cluster <- kmeans(data[, 3:14], centers = n.cluster)
# Cluster ID zum ursprünglichen Datensatz hinzufügen
data$cluster <- fit.cluster$cluster

Nach dem wir unsere Whiskies in 4 Cluster eingeteilt haben, betrachten wir nun deren geschmackliche Eigenschaften. Dazu berechnen wir für jeden Geschmacksindikator innerhalb der Cluster einen Mittelwert.

# Cluster Eigenschaften berechnen
cluster.means <- aggregate(data[, 3:14], by = list(fit.cluster$cluster), FUN = mean)
names(cluster.means)[1] <- "Cluster" 
cluster.means <-melt(cluster.means, id = "Cluster", value.name = "Mean")
names(cluster.means)[2] <- "Taste"

Diese Mittelwerte visuallisieren wir wiederum mit ggplot2:

# Unsere Lieblingsfarben
stat.color <- c("#0A71B4", "#61636B", "#13235B", "#E36929", "#C9D30F", "66B8DC")
# Plot der Geschmacksindikatoren
ggplot(cluster.means, aes(x = Cluster, y = Mean, fill = as.factor(Cluster))) + geom_bar(stat = "identity") + facet_wrap(~ Taste) + scale_x_continuous(breaks = 1:6) + theme_bw() + scale_fill_manual(values = stat.color, name = "Cluster")

plot of chunk unnamed-chunk-11

Wie wir sehen können, empfiehlt sich für Freunde rauchiger Whiskies eindeutig Cluster 2, wohingegen Whiskies aus Cluster 1 einen besonders fruchtig-süßen Charakter besitzen. Liebhaber von Sherry-oder Portcask gelagerten Whiskies werden in Cluster 3 fündig, der sich durch seine wein-ähnlichen Aromen definiert ist. Der “Allrounder”, welcher von allen Geschmacksrichtungen etwas bietet, findet sich in Cluster 4.

# Cluster 1: Süß, floral und fruchtig
cluster.1 <- as.character(droplevels(subset(data, cluster == 1)$Distillery))
# Cluster 2: Rauchig, medizinisch, viel Körper
cluster.2 <- as.character(droplevels(subset(data, cluster == 2)$Distillery))
# Cluster 3: Wein-ähnlich, malzig, nussig und süß
cluster.3 <- as.character(droplevels(subset(data, cluster == 3)$Distillery))
# Cluster 4: Von allem ein bisschen ;)
cluster.4 <- as.character(droplevels(subset(data, cluster == 4)$Distillery))

Empfehlungen anhand der Whisky-Cluster

Wir empfehlen die folgenden Destillerien für Cluster 1 (süß, floral und fruchtig): AnCnoc, ArranIsleOf, Auchentoshan, Aultmore, Benriach, Bladnoch, Bunnahabhain, Cardhu, Craigganmore, Dalwhinnie, Dufftown, GlenElgin, GlenGrant, GlenKeith, GlenMoray, GlenSpey, Glenallachie, Glenfiddich, Glengoyne, Glenkinchie, Glenlossie, Glenmorangie, Inchgower, Linkwood, Loch Lomond, Mannochmore, Miltonduff, RoyalBrackla, Speyburn, Speyside, Strathmill, Tamdhu, Tamnavulin, Teaninich, Tobermory, Tomintoul, Tullibardine

Wir empfehlen die folgenden Destillerien für Cluster 2 (rauchig, medizinisch, viel Körper): Ardbeg, Caol Ila, Clynelish, Lagavulin, Laphroig, Talisker

Wir empfehlen die folgenden Destillerien für Cluster 3 (wein-ähnlich, malzig, nussig und süß): Aberfeldy, Aberlour, Auchroisk, Balmenach, Belvenie, BenNevis, Benrinnes, Benromach, BlairAthol, Craigallechie, Dailuaine, Dalmore, Deanston, Edradour, GlenOrd, Glendronach, Glendullan, Glenfarclas, Glenlivet, Glenrothes, Glenturret, Knochando, Longmorn, Macallan, Mortlach, RoyalLochnagar, Scapa, Strathisla

Wir empfehlen die folgenden Destillerien für Cluster 4 (von allem ein bisschen): Ardmore, Balblair, Bowmore, Bruichladdich, GlenDeveronMacduff, GlenGarioch, GlenScotia, Highland Park, Isle of Jura, Oban, OldFettercairn, OldPulteney, Springbank, Tomatin, Tormore

Whisky-Map of Scotland

Wer sich etwas mit Whiskies auskennt weiß, dass die Region in denen die Destillerien ansässig sind einen entscheidenden Einfluss auf den Charakter des Whiskies haben. So produzieren die Destillerien auf den schottischen Inseln (Islay & Skye) traditionell eher rauchige Whiskies. Im Gegensatz dazu haben Whiskies aus Destillerien an der Nordküste Schottlands (Spreyside) einen eher fruchtig süßen Charakter.

Um diesen Zusammenhang etwas besser zu sehen, lassen wir uns die Destillerien, unterteilt nach Geschmacks-Clustern, auf einer Karte von Schottland anzeigen.

# Als erstes laden wir uns eine Schottlandkarte von stamen.com
map <- get_map(location = c(-4.223918, 57.479072), maptype = "toner", source = "stamen", zoom = 7)
# Diese speichern wir mit ggmap als Objekt mit dem wir später mit ggplot weiter arbeiten können
scotland <- ggmap(map)
# Der von uns genutzte Datensatz beinhaltet bereits die Koordinaten aller Distillen
whiskies.coord <- data.frame(data$Latitude, data$Longitude)
coordinates(whiskies.coord) <- ~data.Latitude + data.Longitude
# Um mit ggmap weiter arbeiten zu können müssen die Koordinaten in Latitue und Longitude übersetzt werden (wir nutzen dazu rgdal und maptools)
proj4string(whiskies.coord) <- CRS("+init=epsg:27700")
whiskies.coord <- spTransform(whiskies.coord, CRS("+init=epsg:4326"))
# Am Ende erstellen wir uns unser Data Frame...
data.map <- data.frame(Cluster = data$cluster,
                       Distillery = data$Distillery,
                       Latitude = whiskies.coord$data.Latitude,
                       Longitude = whiskies.coord$data.Longitude
                       )
# ...u nd fügen die Distillen jeweils mit deren Cluster zur Karte hinzu
scotland + geom_point(data = data.map, aes(x = Latitude, y = Longitude, color = as.factor(Cluster)), size = 6, alpha = .75) + scale_color_manual(values = stat.color, "Cluster") + xlab("Latitude") + ylab("Longitude")

plot of chunk unnamed-chunk-13

Statistischer Zusammenhang zwischen Region und Rauchgehalt

In einem letzten Schritt soll nun die allseits bekannte Hypothese, dass rauchige Whiskies vor allem von den schottischen Inseln kommen, statistisch überprüft werden. Wir nutzen dazu ein sog. lineares Regressionsmodel, das den Effekt einer oder mehrerer Einflussgrößen auf eine Zielgröße schätzt. In unserem Falle ist die Zielgröße des Modells der Rauchgehalt der Whiskies, die Einflussgrößen sind geografische Länge(Longitude) und Breite (Latitude) sowie deren Wechselwirkung. Theoretisch erwarten wir, dass Whiskies im Westen Schottlands (kleinere Werte für Longitude) rauchiger sind als im Osten. Gleichzeitig aber im Norden (größere Werte für Latitude) die Whiskies weniger rauchig sind als im Süden (hier vor allem Isle of Jura und Islay).

# Als erste Verbinden wir die Daten der Geschmacksrichtung mit den Koordinaten in einem Dataframe...
data.map$Smoky <- data$Smoky
# und rechnen dann ein lineares Regressionsmodell
model.1 <- lm(Smoky ~ Longitude * Latitude , data = data.map)
# Als nächstes speichern wir das lm-Objekt aus der Regression in xtable-Objekt damit wir später automatisch eien Tabelle erstellen können
model.1.table <- xtable(model.1, caption = "Table 1: OLS Regression von Längen- und Breitengrad auf Rauchigkeit eines Whiskies")
# Über den "print" Befehl wird die Tabelle dann ausgegeben
print(model.1.table, type = "html")
Table 1: OLS Regression von Längen- und Breitengrad auf Rauchigkeit eines Whiskies
Estimate Std. Error t value Pr(>|t|)
(Intercept) -87.9795 27.5185 -3.20 0.0020
Longitude 1.5430 0.4833 3.19 0.0020
Latitude -18.3527 6.4474 -2.85 0.0056
Longitude:Latitude 0.3159 0.1140 2.77 0.0069


Der geschätze Koeffizient für den Längengrad beträgt 1.54 und ist signifikant (p-value: 0.002). Für den Breitengrad schätzt unser Model einen Wert von -18.35, wiederum signifikant (p-value: 0.006). Je weiter wir uns also in den Osten begeben (zunehmende Latitude), desto weniger rauchig werden die Whiskies im Durchschnitt - das erwarete Ergebnis. Aber, bewegt man sich weiter in den Norden Schottlands (zunehmende Longitude), so steigt die Rauchigkeit der Whiskies statistisch gesehen ebenfalls an. Dieses Resultat wird wahrscheinlich durch die zwei Destillerien aus Cluster 2, die sich sehr weit nördlich befinden beeinflusst (Talisker auf der Isle of Skye sowie Clynelish bei Brora, Highland) beeinflusst.

Für die Interaktion zwischen beiden Prädiktoren schätz unser Model einen Koeffizienten von 0.32 (p-value: 0.007). Da dieser Wert ohne Visualisierung nur schwer interpretierbar ist, visualisieren wir den Koeffizienten in einem Interaktionsplot. Im Folgenden werden zunächst 5 Szenarien definiert, die sich jeweils nur in den Breitengraden, nicht aber den Längengraden unterscheiden. Für jedes Szenario wird dann der Rauchgehalt mit Hilfe unseres Models vorhergesagt.

# Für die Interaktionseffkete erstellen wir zunächst 5 unterschiedliche Szenarien
sc1 <- data.frame(Longitude = seq(55, 59, length.out = 100), Latitude = -6)
sc2 <- data.frame(Longitude = seq(55, 59, length.out = 100), Latitude = -5)
sc3 <- data.frame(Longitude = seq(55, 59, length.out = 100), Latitude = -4)
sc4 <- data.frame(Longitude = seq(55, 59, length.out = 100), Latitude = -3)
sc5 <- data.frame(Longitude = seq(55, 59, length.out = 100), Latitude = -2)
# Mit Hilfe von predict() können wir den Rauchgehalt durch das Modell vorhersagen lassen
y1 <- as.data.frame(predict(model.1, newdata = sc1, interval="prediction")) 
y2 <- as.data.frame(predict(model.1, newdata = sc2, interval="prediction"))
y3 <- as.data.frame(predict(model.1, newdata = sc3, interval="prediction"))
y4 <- as.data.frame(predict(model.1, newdata = sc4, interval="prediction"))
y5 <- as.data.frame(predict(model.1, newdata = sc5, interval="prediction"))
# Schließlich werden die verschiedenen Effekte in Grafiken dargestellt...
p1 <- ggplot(y1, aes(x = seq(55, 59, length.out = 100), y = fit)) + geom_line(color = stat.color[1]) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .3) + ylab("Smoky") + coord_cartesian(ylim = c(0, 5), xlim = c(55, 59)) + xlab("Longitude") + theme_bw() + ggtitle("Latitude = -6")
p2 <- ggplot(y2, aes(x = seq(55, 59, length.out = 100), y = fit)) + geom_line(color = stat.color[1]) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .3) + ylab("Smoky") + coord_cartesian(ylim = c(0, 5), xlim = c(55, 59)) + xlab("Longitude") + theme_bw() + ggtitle("Latitude = -5")
p3 <- ggplot(y3, aes(x = seq(55, 59, length.out = 100), y = fit)) + geom_line(color = stat.color[1]) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .3) + ylab("Smoky") + coord_cartesian(ylim = c(0, 5), xlim = c(55, 59)) + xlab("Longitude") + theme_bw() + ggtitle("Latitude = -4")
p4 <- ggplot(y4, aes(x = seq(55, 59, length.out = 100), y = fit)) + geom_line(color = stat.color[1]) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .3) + ylab("Smoky") + coord_cartesian(ylim = c(0, 5), xlim = c(55, 59)) + xlab("Longitude") + theme_bw() + ggtitle("Latitude = -3")
p5 <- ggplot(y5, aes(x = seq(55, 59, length.out = 100), y = fit)) + geom_line(color = stat.color[1]) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .3) + ylab("Smoky") + coord_cartesian(ylim = c(0, 5), xlim = c(55, 59)) + xlab("Longitude") + theme_bw() + ggtitle("Latitude = -2")
# und ausgegeben
grid.arrange(p1, p2, p3, p4, p5, ncol= 2, nrow = 3)

plot of chunk unnamed-chunk-17

Zusammenfassung und Ausblick

Bereits heute können durch die Integration von knitr in RStudio vielfälltige Reporting-Aufgaben automatisiert werden. In einer Zeit in der sich Daten innerhalb von Bruchteilen von Sekunden ändern, ist dies ein notwendiger Schritt um effizient Analysen durchzuführen und Ergebnisse zu publizieren.

In diesem Beitrag haben wir einige der Möglichkeiten, die knitr und RStudio bieten, gezeigt. Unter anderem können Variablen durch Verweise im Fließtext dynamisch eingebunden werden. Gleiches gilt auch für Tabellen, die sich zum Beispiel mit Hilfe des “xtable-Packages” automatisch erzeugen und ausgeben lassen. Schlussendlich bietet “ggplot2” dem Benutzer die Möglichkeit state-of-the-Art Grafiken zu produzieren und diese direkt in den Text zu integrieren.

Damit sind aber noch lange nicht alle Möglichkeiten die das automatisierte Reporting in R bietet ausgenutzt. So lassen sich zum Beispiel mit Hilfe von shiny dynamische Applikationen erstellen und direkt in Reports integrieren. Schlussendlich kann mit Hilfe von latex vollständig benutzerspezifischer Output generiert werden.

In eigener Sache

Als Whisky-Liebhaber hat uns die Analyse der Daten natürlich besonders viel Spaß gemacht. Interessierte, die über genauerer Daten, z.B. über einzelne Whiskies und nicht nur Destillerien, verfügen, dürfen Sie gerne unter +49-(0)69-6783-0675-1 oder unter info@statworx.de melden. Wir freuen uns auf euer Feedback!

Folgen Sie uns auch auf: