Einleitung

Quantile und die daraus abgeleiteten Box-Plots als grafische Zusammenfassung sind geeignete Mittel, um die Streuung der Daten zu charakterisieren. In diesem Tutorial wird gezeigt, wie die 5-Statistiken eines Box-Plots in R numerisch bestimmt werden. Die hier verwendeten Daten stammen aus dem mtcars-Datensatz, welcher im dataset-Package zu finden ist.

Zu Beginn werden die Libraries geladen. Mit ggplot2 können komplexe Grafiken sehr leicht erzeugt werden. Die Funktion basiert auf der “Grammar of Graphics”. Die Library xtable wird ebenfalls benötigt.

sessionInfo()
## R version 3.2.0 (2015-04-16)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
## 
## locale:
## [1] de_CH.UTF-8/de_CH.UTF-8/de_CH.UTF-8/C/de_CH.UTF-8/de_CH.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] magrittr_1.5    formatR_1.2     tools_3.2.0     htmltools_0.2.6
##  [5] yaml_2.1.13     stringi_0.4-1   rmarkdown_0.6.1 knitr_1.10.5   
##  [9] stringr_1.0.0   digest_0.6.8    evaluate_0.7
library(ggplot2)
library(xtable)

Datensatz mtcars

Der mtcars-Datensatz beinhaltet Informationen zu bestimmten Wagenmodellen. In diesem Tutorial legen wir den Fokus auf die Anzahl gefahrenen Meilen pro Gallone der verschiedenen Wagenmodellen. Deren Motorengrösse ist mit der Anzahl Zylinder bestimmt. Die Anzahl Meilen pro Gallone ist in der Variable mpg erfasst. Die Anzahl Zylinder ist in der Variable cyl enthalten. Die übrigen Variablen im Datensatz werden in diesem Tutorial nicht verwendet.

head(mtcars)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Die Variable cyl wird in eine factor-Variable konvertiert. Zusätzlich wird eine neue Variable labels kreiert. Die letzt genannte Variable speichert die Namen der verschiedenen Wagenmodelle.

mtcars$cyl <- factor(mtcars$cyl)
mtcars$labels <- row.names(mtcars)

Eine Zusammenfassung des Datensatzes gibt die Gelegenheit sich mit den verschiedenen Variablen vertraut zu machen.

summary(mtcars)
##       mpg        cyl         disp             hp             drat      
##  Min.   :10.40   4:11   Min.   : 71.1   Min.   : 52.0   Min.   :2.760  
##  1st Qu.:15.43   6: 7   1st Qu.:120.8   1st Qu.: 96.5   1st Qu.:3.080  
##  Median :19.20   8:14   Median :196.3   Median :123.0   Median :3.695  
##  Mean   :20.09          Mean   :230.7   Mean   :146.7   Mean   :3.597  
##  3rd Qu.:22.80          3rd Qu.:326.0   3rd Qu.:180.0   3rd Qu.:3.920  
##  Max.   :33.90          Max.   :472.0   Max.   :335.0   Max.   :4.930  
##        wt             qsec             vs               am        
##  Min.   :1.513   Min.   :14.50   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :3.325   Median :17.71   Median :0.0000   Median :0.0000  
##  Mean   :3.217   Mean   :17.85   Mean   :0.4375   Mean   :0.4062  
##  3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :5.424   Max.   :22.90   Max.   :1.0000   Max.   :1.0000  
##       gear            carb          labels         
##  Min.   :3.000   Min.   :1.000   Length:32         
##  1st Qu.:3.000   1st Qu.:2.000   Class :character  
##  Median :4.000   Median :2.000   Mode  :character  
##  Mean   :3.688   Mean   :2.812                     
##  3rd Qu.:4.000   3rd Qu.:4.000                     
##  Max.   :5.000   Max.   :8.000

Box-Plot

Mit der ggplot() Funktion wird ein Box-Plot erstellt. Die Grafik zeigt die Anzahl gefahrenen Meilen pro Gallone der verschieden grosse Motoren an. Die Motorengrösse ist mit der Anzahl Zylinder angegeben. Insgesamt sind drei Motorengruppen mit 4, 6 und 8 Zylindern im Datensatz vertreten.

Der Box-Plot zeigt, dass Motoren mit einer grösserer Anzahl Zylindern weniger Meilen pro Gallone fahren. D.h. je grösser der Motor, desto weniger Meilen pro Gallone werden gefahren. Was intuitive Sinn macht, da ein grösserer Motor mehr Benzin verbraucht als ein kleinerer Motor.

p <- ggplot(data = mtcars, aes(x = cyl, y = mpg, fill = cyl))
p <- p + geom_boxplot() + 
  ggtitle("Car Milage Data") + 
  labs(x = "Number of Cylinders", y = "Miles Per Gallon") + 
  scale_fill_discrete(name = "Cylinders")
p

Die Gruppe mit einer Motorengrösse von 8 Zylinder zeigt zwei Datenpunkte, die ausserhalb der blauen Box zu liegen kommen. Diese Datenpunkte sind als potenzielle Ausreisser identifiziert worden. Sind es wirklich nur zwei Datenpunkte? Hierfür wird die Funktion add.outlier() verwendet. Mit der Funktion werden die identifizierten Ausreisser in der Grafik mit den jeweiligen Werten oder Labels gekennzeichnet.

add.outlier <- function(p,labvar = as.character(p$mapping$y)){
  df <- data.frame(y = with(p$data,eval(p$mapping$y)),
                   x = with(p$data,eval(p$mapping$x)))
  
  df.l <- split(df,df$x)
  
  mm <- Reduce(rbind, lapply(df.l,FUN = function(df){
    data.frame(y = df$y[df$y <= (quantile(df$y)[2] - 1.5 * IQR(df$y)) | df$y >= (quantile(df$y)[4] + 1.5 * IQR(df$y))],
               x = df$x[df$y <= (quantile(df$y)[2] - 1.5 * IQR(df$y)) | df$y >= (quantile(df$y)[4] + 1.5 * IQR(df$y))]
    )})
  )
  
  
  mm$x <- factor(mm$x,levels=sort(as.numeric(as.character(unique(p$data[,as.character(p$mapping$x)])))),
                 labels = levels(p$data[,as.character(p$mapping$x)])
  )
  
  names(mm) <- c(as.character(p$mapping$y),as.character(p$mapping$x))
  mm <- merge(p$data[,c(names(mm),labvar)],mm)
  
  p + geom_text(data=mm,
                aes_string(label=labvar),
                vjust = -0.5)
}

Hier der Box-Plot mit den Ausreisser-Werten.

add.outlier(p)

Hier der Box-Plot mit den Ausreisser-Labels. Dabei ist ersichtlich, dass der Box-Plot der 8 Zylinder Motoren unterhalb der Box gleich zwei potenzielle Ausreisser identifiziert hat. Die Ausreisser überlappen sich. Es handel sich um die zwei Wagenmodelle “Cadillac Fleetwood” und “Lincoln Continental”.

add.outlier(p,"labels")

row.names(mtcars)[mtcars$cyl == 8 & mtcars$mpg == 10.4]
## [1] "Cadillac Fleetwood"  "Lincoln Continental"

Interquartilsabstand IQR

Eine Masszahl für die Streuung der Daten ist der Interquartilsabstand (\(IQR\) := “interquartile range”). Der IQR berechnet sich aus der Differenz des 3. und 1. Quartils (\(IQR = x_{0.75} - x_{0.25}\)). In der Grafik entspricht der IQR der Distanz der Box.

Der IQR ist resistent gegen Ausreisser, da die Quartile nicht von der Lage der Daten links von \(x_{0.25}\) und rechts von \(x_{0.75}\) beeinflusst werden. Die Quartile der Motoren mit 8 Zylindern werden mit der Funktion quantile() bestimmt. Hier die empirisch geschätzten Quantile:

quantile(mtcars$mpg[mtcars$cyl == 8])
##    0%   25%   50%   75%  100% 
## 10.40 14.40 15.20 16.25 19.20

Ein Quartil ist ein bestimmtes Quantil. Das 25%-Quantil und das 75%-Quantil entsprechen dem 1. bzw. 3. Quartil. Der Median \(\hat{x}_{0.5}\) ist 15.2. Die Quartile \(\hat{x}_{0.25}\) und \(\hat{x}_{0.75}\) entsprechen 14.4 und 16.25.

Whisker

Die Linien ausserhalb der Box sind die sogenannten Whisker. Die Whisker Linien sind keine Vertrauensintervalle und sie entsprechen auch nicht dem Minimum und/oder Maximum der Daten. Die Whisker werden immer mit Hilfe einer Faustregel bestimmt. Die Whiskergrenzen laden immer auf einen Datenpunkt im Datensatz. Alle Datenpunkte ausserhalb des Whiskers werden als potenzielle Ausreisser-Kandidaten identifiziert. Um die Whisker zu bestimmen, muss die Unter- und Obergrenze definiert werden. Dafür wird eine Faustregel verwendet. Basierend auf der Unter- und Obergrenze werden auch die Ausreisser bestimmt.

Eine Faustregel zur Identifikation von potentiellen Ausreissern ist: Bilde den inneren “Zaun” mit der Untergrenze \(z_{u} = x_{0.25} - 1.5 * IQR\) und der Obergrenze \(z_{o} = x_{0.75} + 1.5 * IQR\). Daten kleiner als \(z_{u}\) und grösser als \(z_{o}\) sind dann Ausreisser-Kandidaten, die genauer zu inspizieren sind. Mit Hilfe der oben berechneten Quartile kann die Unter- und Obergrenze exakt bestimmt werden. Nach unserer eingeführten Faustregel lauten diese folgendermassen: \(\hat{IQR} = (16.25 - 14.4)\), \(\hat{z}_{u} = 14.4 - 1.5 * \hat{IQR} = 11.625\) und \(\hat{z}_{o} = 16.25 + 1.5 * \hat{IQR} = 19.025\).

Eine Alternative um direkt die Whiskergrenzen \(w_{u}\) und \(w_{o}\) der Motoren mit 8 Zylindern zu bestimmen, ist die Funktion boxplot() zu verwenden.

b <- boxplot(mpg ~ cyl, data = mtcars[mtcars$cyl == 8,])

b$stats
##      [,1] [,2] [,3]
## [1,]   NA   NA 13.3
## [2,]   NA   NA 14.3
## [3,]   NA   NA 15.2
## [4,]   NA   NA 16.4
## [5,]   NA   NA 19.2

Dabei ist der Median \(\hat{x_{box}}_{0.5}\) 15.2. Die Quartile \(\hat{x_{box}}_{0.25}\) und \(\hat{x_{box}}_{0.75}\) entsprechen 14.3 und 16.4. Die Whiskergrenzen \(\hat{w_{box}}_{u}\) und \(\hat{w_{box}}_{o}\) sind 13.3 und 19.2.

Im Unterschied zur ggplot() Funktion wird hier der Datenpunkt 19.2 nicht als Ausreisser identifiziert. Anscheinend verwendet die boxplot() Funktion eine andere Faustregel zur Ermittlung der Unter- und Obergrenze. Wie dieser Algorithmus genau funktioniert, ist zu ergründen. Was aber nicht Bestanteil des Tutorial ist. Die Faustregel für die boxplot() Funktion ist aber sehr ähnlich wie die bereits eingeführte Faustregel. Anstatt die empirisch geschätzten Quartile werden die beobachteten Daten verwendet, die nahe an den empirisch geschätzten Quartilen liegen:

\(\hat{IQR} = (16.4 - 14.3)\), \(\hat{z}_{u} = 14.3 - 1.5 * \hat{IQR} = 11.15\) und \(\hat{z}_{o} = 16.4 + 1.5 * \hat{IQR} = 19.55\)

Die Unter- und Obergrenzen wurden einerseits mit der eingeführten Faustregel mittels der empirisch geschätzen Quartile und andererseits mit Hilfe der beobachteten Daten ermittelt. Wie aus den Berechnungen hervorgeht, sind die Unter- und Obergrenzen nicht in beiden Fällen gleich. Im nächsten Kapitel wird der Fokus nochmals auf diesen Unterschied gelegt. Dabei wird erklärt, weshalb die beiden Varianten (ggplot() und boxplot()) unterschiedliche Ausreisser-Kandidaten im Datensatz identifizieren.

Die 5-Statistiken

In der folgenden Tabelle werden die ermittelten Unter- und Obergrenzen mit der ggplot() Funktion und der boxplot() Funktion gezeigt.

ggplot boxplot
Obere G. 19.02 19.55
3.Q 16.25 16.40
Median 15.20 15.20
1.Q. 14.40 14.30
Untere G. 11.62 11.15

Wenn man die Obergrenze genauer unter die Lupe nimmt, ist ersichtlich, dass sie im Fall der boxplot() Funktion höher ist als in der ggplot() Funktion. Dies erklärt weshalb die boxplot() Funktion den Datenpunkt 19.2 nicht als potenziellen Ausreisser-Kandidaten identifiziert. Um nochmals kurz in Erinnerung zu rufen, Datenpunkte, die ausserhalb der Unter- oder Obergrenze liegen, werden als Ausreisser identifiziert. Der Datenpunkt 19.2 liegt nicht ausserhalb der Obergrenze bei 19.55. Der Datenpunkt 19.2 wird also nicht als Ausreisser identifiziert. Im Fall der ggplot() Funktion, wo die Obergrenze bei 19.025 liegt, ist der Datenpunkt 19.2 ausserhalb der Obergrenze, so dass er als potenziellen Ausreisser identifiziert wird.

Zur Veranschaulichung der beiden unterschiedlichen “Faustregeln” wird der sortierte Datensatz der Motoren mit 8 Zylinder gezeigt. Zusätzlich werden die Box-Plots mit deren horizontalen Unter- und Obergrenzen und deren 1. und 3. Quartilslinien dargestellt. Die Whsikergrenzen spannen die Streuung der Daten auf. Vom Median aus gesehen sind die Whiskergrenzen im Allgemeinen höchstens so lang wie die Unter- und Obergrenzen.

Der Datensatz der Motoren mit 8 Zylindern.

##  [1] 10.4 10.4 13.3 14.3 14.7 15.0 15.2 15.2 15.5 15.8 16.4 17.3 18.7 19.2

Der ggplot() Box-Plot mit den berechneten Unter- und Obergrenzen. Die Quartile sind empirisch geschätzt.

Der boxplot() Box-Plot mit den berechneten Unter- und Obergrenzen. Die Quartile basierend auf die beobachteten Daten.