Gerade mit Referenz-Steigung in einen doppelt logarithmischen Plot in R eintragen

Ich stelle viele meiner Daten in Plots mit zwei logarithmischen Achsen dar. Oftmals kommt es in solchen Plots auf die Steigung der Kurve an, also den Exponenten eines Power laws. In solchen Fällen zeichne ich, zusätzlich zu den Daten, eine oder mehrere Geraden mit einer Referenz-Steigung ein. Dadurch kann der Betrachter abschätzen, in welchen Bereichen die Daten eher zu der einen oder der anderen Steigung tendieren. Ich fand es allerdings immer sehr mühsam diese Steigungen einzuzeichnen. Grund: Ein Teil der Anfangs- oder Endkoordinaten musste im logarithmischen Raum berechnet werden. Daher habe ich mir eine Funktion geschrieben, die mir anhand einer Anfangs-Koordinate, der Steigung, sowie einem Teil der Endkoordinate, den fehlenden Wert berechnet und die Referenz-Gerade in den Plot einzeichnet.

  1. Zunächst generiere ich Daten, die im doppelt-logarithmischen Plot eine Gerade ergeben:
    x = 1:1000
    y = x^-1.5
  2. Diese werden dann mit schön formatierten Achsen geplottet, wie im Artikel „Logarithmische Achsen in R-Plots formatieren“ beschrieben.
    library(package = "sfsmisc")
    #png("log-gerade.png", width=480, height=360)
    #par(cex =1.5, mar=c(4,4,0,0)+.01)
     
    plot(x, y, log="xy", xaxt="n", yaxt="n")
    sfsmisc::eaxis(side=1, at=10^c(0:3))
    sfsmisc::eaxis(side=2, at=10^c(-4,-2,0))
     
    #dev.off()
    Darstellung einer Potenz-Funktion in einem Plot mit zwei logarithmischen Achsen.
    Darstellung einer Potenz-Funktion in einem Plot mit zwei logarithmischen Achsen.
  3. Der Exponent der von mir verwendeten Potenzfunktion ist -1.5. Daher ist auch die Steigung der Geraden m = -1.5. Nun möchte ich zusätzlich eine Gerade mit der Steigung m = -1 in den Plot eintragen. Diese soll bei am Punkt (x1 = 101, y1 = 10-1) beginnen. Enden soll die Gerade bei x2 = 2*102. Gesucht ist also die Y-Koordinate des Endpunktes (y2).

    Nun kann man:

    1. Alle vorhandenen x- und y-Werte logarithmieren (da beide Achsen logarithmisch sind).
    2. Diese Werte können in die Gleichung für die Berechnung der Steigung eingesetzt werden.
    3. Nun muss die Gleichung nach x2 aufgelöst werden (oder: die Gleichung auflösen lassen).
    4. Im Anschluss muss der x2-Wert in den Exponenten genommen werden.
    5. Über den Plot-Befehl kann nun eine Gerade zwischen den beiden Punkten (x1, y1) und (x2, y2) eingezeichnet werden.

    Diese Vorgehensweise habe ich in der Funktion refline() zusammengefasst:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    
    refline <- function(x1=NA, y1=NA, x2=NA, y2=NA, m=1, add=T, ...) {
      # Welchen Achsen (des letzten plots) sind logarithmisch?
      logX = par("xlog")
      logY = par("ylog")
     
      # Werte der logarithmischen Achsen umrechnen
      if(!is.na(x1) && logX) {x1 = log(x1)}
      if(!is.na(x2) && logX) {x2 = log(x2)}
      if(!is.na(y1) && logY) {y1 = log(y1)}
      if(!is.na(y2) && logY) {y2 = log(y2)}
     
      # Bestimmung der fehlenden Variablen
      if(is.na(x1)) {
        x1 = (y1 - y2 + m * x2) / m
        print("X1 wird berechnet")
      } else if(is.na(x2)) {
        x2 = (y2 - y1 + m * x1) / m
        print("X2 wird berechnet")
      } else if(is.na(y1)) {
        y1 = y2 + m * (x1 - x2)
        print("Y1 wird berechnet")
      } else if(is.na(y2)) {
        y2 = y1 + m * (x2 - x1)
        print(paste("Y2 wird berechnet", y2))
      }
     
      # Werte der logarithmischen Achsen umrechnen
      if(!is.na(x1) && logX) {x1 = exp(x1)}
      if(!is.na(x2) && logX) {x2 = exp(x2)}
      if(!is.na(y1) && logY) {y1 = exp(y1)}
      if(!is.na(y2) && logY) {y2 = exp(y2)}
     
      # Einzeichnen
      if(add) {
        lines(x=c(x1, x2), y=c(y1, y2), ...)
      }
     
      # Ausgabe
      data.frame(
        x1, y1, x2, y2, m
        , xmin=min(x1, x2), xmax=max(x1, x2)
        , ymin=min(y1, y2), ymax=max(y1, y2)
      )
    }

    Die Besonderheiten der Funktion refline() sind:

    • Es wird automatisch die fehlende Koordinate bestimmt.
    • Es wird erkannt, welche der Achsen logarithmisch sind.
    • Wichtige Werte der Geraden werden bei Aufruf zurückgegeben. Im Anschluss können sie verwendet werden, z.B. um Texte an die Enden der Geraden zu schreiben.
    • Die Gerade wird (auf Wunsch: add=T) in den vorigen Plot eingetragen.
    • Grafik-Parameter werden an den Plot-Befehl weitergeleitet (...).

    Nun kann refline() genutzt werden, um die oben genannte Gerade (sowie eine weitere Gerade) in den Plot einzuzeichnen. Die Variablen rl1 und rl2 dienen dazu, die Anfangs- und End-Punkte zu speichern. Mit den gespeicherten Werten werden Anschluss die Texte positioniert.

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    
    #png("log-gerade_2.png", width=480, height=360)
    #par(cex =1.5, mar=c(4,4,0,0)+.01)
     
    plot(x, y, log="xy", xaxt="n", yaxt="n")
    sfsmisc::eaxis(side=1, at=10^c(0:3))
    sfsmisc::eaxis(side=2, at=10^c(-4,-2,0))
     
    rl1 = refline(x1=1e1, y1=1e-1, x2=2e2, m=-1, col="red")
    text(x=rl1$x2, y=rl1$y2, labels=paste0(" Slope: ", rl1$m), adj=0, col="red")
     
    rl2 = refline(x1=8e0, y1=1e-2, x2=1e2, m=-2, col="blue")
    text(x=rl2$x2, y=rl2$y2, labels=paste0(" Slope: ", rl2$m), adj=0, col="blue")
     
    #dev.off()
    Darstellung einer Potenzfunktion mit einer Steigung von -1.5. Zusätzlich wurden Referenz-Steigungen von -1 und -2 eingezeichnet.
    Darstellung einer Potenzfunktion mit einer Steigung von -1.5. Zusätzlich wurden Referenz-Steigungen von -1 und -2 eingezeichnet.
  4. Hinweis: Als Parameter für die Funktion refline() muss ein vollständiges und ein unvollständiges Koordinaten-Tupel, sowie die gewünschte Steigung angegeben werden. Der Punkt (x2, y2) muss dabei nicht immer rechts vom Punkt (x1, y1) liegen. Vielmehr hängt es von den gewählten Werten ab, wie die Punkte zueinander liegen. Die Werte x1 und y1 bestimmen gemeinsam immer den einen Punkt und die Werte x2 und y2 gemeinsam den anderen Punkt.
  5. Wie ich oben bereits erwähnt habe, funktioniert die Berechnung des fehlenden Punktes der Gerade auch in Plots mit zwei linearen, oder nur einer logarithmischen Achse:
    Im Falle z.B. einer logarirhmischen Y-Achse werden dann nur die Y-Koordinaten logarithmiert. Hier einige Beispiele:

    # Linear
    plot(x, y)
    r1 = refline(x1=0, x2=500, y1=1, m=-0.0005, col="blue")
    text(x=r1$xmax, y=r1$ymin, labels=paste0(" Steigung: ", r1$m), adj=0)
     
    # Y-logarithmisch
    plot(x, y, log="y", yaxt="n")
    sfsmisc::eaxis(side=2, at=10^c(-4,-2,0))
    r2 = refline(x1=10, x2=500, y1=1, m=-0.002, col="red")
    text(x=r2$xmax, y=r2$ymin, labels=paste0(" Steigung: ", r2$m), adj=0)

Zeichenketten in R zusammenfügen

Das Zusammenfügen von Zeichenketten (Konkatenieren von Strings) ist eine Standard-Aufgabe in allen Programmiersprachen. Auch bei der Programmierung in R wird die Funktion benötigt, wenn man z.B. Text dynamisch erzeugen lassen, oder den Wert einer Variablen an einen Text anfügen möchte. Die Funktion, die diese Aufgabe in R übernimmt, ist meiner Meinung nach etwas außergewöhnlich benannt. Mit diesem Artikel möchte ich mich selbst daran erinnern, dass sie paste() heißt.

Hier ein paar Beispiele, wie Zeichenketten (Strings) in verschiedenen Programmiersprachen zusammengefügt werden:

R
Die Funktion, mit der in R Strings zusammengefügt werden, heißt paste() (vom englischen „paste“ für „zusammenkleben“). Beispiel:

daten = c(1:10)
plot(daten, main=paste("Plot von",length(daten),"Werten."))

Mit der Funktion paste() werden durch Kommata getrennte Werte zusammengefügt. Dabei wird jedesmal ein Leerzeichen eingefügt. Wenn man die Zeichenketten ohne Leerzeichen zusammenfügen möchte, kann man entweder den Parameter sep="" setzen, oder die Funktion paste0() verwenden:

paste(sep="", "eins", "zwei") == paste0("eins", "zwei")
Java
In Java können Zeichenketten mit dem +-Operator verkettet werden:

String var = "eins";
var = var + "zwei" + "drei";
System.out.println(var);

oder mit dem +=-Operator:

String text = "Anfang";
text += "Ende";
System.out.println(text);

Sollte man sehr häufig Strings zusammenfügen, so empfiehlt sich die Verwendung eines StringBuffer-Objektes mit der Methode append() (siehe auch ChuckAndWayne-Blog):

StringBuffer b = new StringBuffer();
for (int i = 0; i<1000000; i++) {
	b.append("langer Text\n");
}
System.out.println(b);
bash
In der Bash können Variablen und Strings einfach direkt nacheinander geschrieben werden, um sie zusammenzufügen.

a="A"
b=$a"b"
echo $b
 
c="${b}c"
echo $c
 
echo $c "Bei echo können" "Leerzeichen zwischen" "den Strings sein."
 
d=$a$b
echo $d

Die bash unterstützt auch den +=-Operator:

c="Vorne"
c+=Hinten
echo $c

Für weitere Details siehe StackOverflow.

Um ganze Dateien zusammenzufügen, kann man den Befehl cat verwenden:

echo -e "Zeile 1, Datei 1\nZeile 2, Datei 1" > datei1.txt
echo -e "Zeile 1, Datei 2\nZeile 2, Datei 2" > datei2.txt
cat datei1.txt datei2.txt
PHP
In PHP können Strings mit Hilfe von . oder .= zusammengefügt werden:

<?php
$text = "Mein Text ";
$text .= "wird immer" . " " . "länger!";
echo $text;
?>
JavaScript
In JavaScript können Strings, ähnlich wie in Java, mit Hilfe des + und des +=-Operators aneinander gehängt werden:

<script>
var msg = "Nachricht";
msg += " wird";
msg = msg + " immer" + " länger";
alert(msg + "!");
</script>
C#
In C# können Strings mit Hilfe der Funktionen string.concat and string.join zusammengefügt werden. Siehe dazu den Blog von Alex James Brown.

Linientypen, Symbole und Farben in R-Plots

R bietet verschiedene Möglichkeiten, seine Datenreihen in Plots unterscheidbar zu machen. Dazu gibt man in den Plot-Befehlen den zu ändernden Parameter, sowie einen Zahlenwert an. Ich vergesse allerdings immer, welche Zahl zu welcher Farbe oder zu welchem Linientyp gehört. Daher liste ich in diesem Artikel Zahlenwerte für folgende Parameter auf:

  1. Typen von Linien (lty).
    Beispiel: plot(c(1:10), type="o", lty=3).
  2. Farben (col).
    Beispiel: plot(c(1:10), type="o", col=4).
  3. Typen von Punkten / Sybole (pch; point character).
    Beispiel: plot(c(1:10), type="o", pch=5).

Linientypen

In R gibt es sechs verschiedene Typen von Linien. Die Zahlen für die Linientypen werden mit einer Periode von sechs „recycled“. Dadurch entspricht der Typ sieben (7) wieder dem Typen eins (1), eine Linie vom Typ acht (8) sieht aus wie eine Linie vom Typen zwei (2) und so weiter. Für jeden der Linien-Typen gibt es auch einen Namen (siehe auch Beschreibung des Parameters lty in der R-Hilfe zu par(): ?par):

  • 0: „blank“; unsichtbare Linie (=> wird nicht gezeichnet)
  • 1: „solid“
  • 2: „dashed“
  • 3: „dotted“
  • 4: „dotdash“
  • 5: „longdash“
  • 6: „twodash“
In R gibt es sechs verschiedene Typen von Linien, hier abwechselnd in rot und grün dargestellt. Jeder der Linientypen ist in 3 verschiedenen Größen (lwd=1-3) dargestellt. Die Werte für die Linientypen wiederholen sich in einer Periode von sechs. Daher entspricht die Angabe lty=7 der Angabe lty=1, die Angabe lty=8 entspricht der Angabe lty=2, und so weiter.
In R gibt es sechs verschiedene Typen von Linien, hier abwechselnd in rot und grün dargestellt. Jeder der Linientypen ist in 3 verschiedenen Größen (lwd=1-3) dargestellt. Die Werte für die Linientypen wiederholen sich in einer Periode von sechs. Daher entspricht die Angabe lty=7 der Angabe lty=1, die Angabe lty=8 entspricht der Angabe lty=2, und so weiter.

Das oben stehende Bild habe ich mit folgendem Code erzeugt:

1
2
3
4
5
6
7
8
# 6 Linientypen
png("linientypen.png")
par(mar=c(0,4,0,0)+.1, cex=1.5)
plot(1,1, type="n", ylim=c(0, 7), xaxt="n", xlab="", ylab="Linentyp (lty)", las=1)
abline(h=c(0:7),     lty=c(0:7), lwd=3, col=c(2,3))
abline(h=c(0:7)-0.2, lty=c(0:7), lwd=2, col=c(2,3))
abline(h=c(0:7)-0.4, lty=c(0:7), lwd=1, col=c(2,3))
dev.off()

Farben

In der Standard-Palette von R befinden sich acht verschiedene Farben. Dazu kommt noch der Wert 0, bei dem die Hintergrundfarbe als Farbe verwendet wird. Die Farben werden mit einer Periode von acht „recycled“. Durch diese fortlaufende Wiederholung, entspricht die Farbe neun (9) der Farbe eins (1), die Farbe zehn (10) wieder der Farbe zwei (2) und so weiter. Anstatt einer Zahl kann man auch direkt den RGB-Farbwert angeben (z.B. plot(1, 1, col="#000000") anstatt plot(1, 1, col=1) ). Weitere Details gibt es im Bereich „Color Specification“ in der R-Hilfe zu den Grafik-Parametern (?par).

In R gibt es acht verschiedene Farben. Hinzu kommt der Wert 0, bei dem die Hintergrundfarbe als Farbe verwendet wird (siehe diagonale Linie). Mit einer Periodenlänge von 8, werden die Farben wiederholt (wiederholte Farben sind mit einer dünnen Linie dargestellt).
In R gibt es acht verschiedene Farben. Hinzu kommt der Wert 0, bei dem die Hintergrundfarbe (hier: #eeeeee) als Farbe verwendet wird (siehe diagonale Linie). Mit einer Periodenlänge von 8, werden die Farben wiederholt (wiederholte Farben sind mit einer dünnen Linie dargestellt).

Für das vorangegangene Bild habe ich folgenden Code verwendet:

1
2
3
4
5
6
7
# 8 LinienFarben
png("linienfarben.png")
par(mar=c(0,4,0,0)+.1, cex=1.5, bg="#eeeeee")
plot(1,1, type="n", ylim=c(0, 11), xlim=c(0, 11), xaxt="n", xlab="", ylab="Linenfarbe (col)", las=1, lab=c(1,11,0))
abline(h=c(0:11), lty=1, lwd=c(rep(10,9), rep(1,9)), col=c(0:11))
abline(a=0, b=1, lwd=10, col=0)
dev.off()

Punkttypen / Symbole

Um den Punkttyp zu definieren, kann entweder ein einzelnes Zeichen oder eine Zahl angeben werden. In der Hilfe zum Befehl par() (?par oder ?graphics::par) steht dazu:

Note that only integers and single-character strings can be set as a graphics parameter (and not NA nor NULL)

Eine detaillierte Liste gibt es in der R-Hilfe zur Funktion points() (?points; Bereich „‚pch‘ values“). Hier eine kurze Zusammenfassung:

  • Zeichen 018: S-kompatible Vektor Symbole
  • Zeichen 1925: Weitere R Vektor Symbole. Die Zeichen 2125 können mit Hilfe von bg innerhalb des plot-Befehls (hier: points()) eingefärbt werden. Achtung: Die Angabe von bg innerhalb von points() hat hier eine andere Auswirkung als die Angabe von bg innerhalb von par() (z.B. par(bg="#ff00ff")). Durch letzteres setzt man die Hintergrundfarbe der gesamten Plot-Fläche.
  • Zeichen 2631: Werden ignoriert
  • Zeichen 32127: ASCII Zeichen

In der folgenden Grafik werden die die ersten 26 Zeichen dargestellt (+ zwei negative Werte).

Punkt-Typen für "pch"-Werte von -2 bis +25. Der Innernraum der Zeichen 21 bis 25 wurde durch die Angabe " bg='red' " innerhalb des Befehls "points()" rot eingefärbt.
Punkt-Typen für „pch“-Werte von -2 bis +25. Der Innernraum der Zeichen 21 bis 25 wurde durch die Angabe “ bg=’red‘ “ innerhalb des Befehls „points()“ rot eingefärbt.

Das oben stehende Bild wurde mit folgenden Befehlen erzeugt:

1
2
3
4
5
6
7
8
9
# 25 Punkt-Typen
png("punkttypen.png")
par(mar=c(0,4,0,0)+.25, cex=1.5)
plot(1, 1, type="n", ylim=c(-2, 25), xlim=c(0, 10), xaxt="n", xlab="", ylab="Punkttyp (pch)", las=1, lab=c(2,15,0))
abline(h=c(-2:25), lty=3, lwd=1, col="#cccccc")
for (typ in -2:25) {
  points(x=c(0:10), y=c(rep(typ,11)), pch=typ, bg="red")
}
dev.off()

Logarithmische Achsen in R-Plots formatieren

Im Gegensatz zu gnuplot, werden in R-Plots die logarithmische Achsen nicht sehr schön formatiert. Das Paket SFSmisc schafft hier Abhilfe. In diesem Artikel beschreibe ich Schritt für Schritt, wie man ansehnliche logarithmische Achsen für R-Plots erstellt.

  1. Als Beispiel sollen Daten geplottet werden, die in einem doppelt logarithmischen Plot eine Gerade ergeben:
    x = 1:1000
    y = x^-1.5
  2. Nun werden diese Daten mit dem Befehl plot() in einem Koordinatensystem mit logarithmischer X- und Y-Achse (log="xy") dargestellt:
    plot(x, y, log="xy")

    Dabei kann man erkennen, dass es keine kleinen Striche zur Unterteilung der Bereiche zwischen den einzelnen Größenordnungen gibt:

    Plot mit logarithmischen Achsen. Standardmäßig werden keine Striche zwischen den Größenordnungen eingezeichnet.
    Plot mit logarithmischen Achsen. Standardmäßig werden keine Striche zwischen den Größenordnungen eingezeichnet.
  3. Es ist aber durchaus üblich, den Bereich zwischen zwei Größenordnungen mit 8 kleinen Strichen zu unterteilen. Diese Striche zeigen an, an welchen Stellen 20%, 30%, 40%, 50%, 60%, 70%, 80% oder 90% der folgenden Größenordnung erreicht wurden. Um das in R zu bewerkstelligen, muss zunächst das Paket SFSmisc installiert werden (siehe dazu auch: Funktion aus einem bestimmten R-Paket laden):
    # Installation
    install.packages("sfsmisc", dependencies=T)
    # Paket laden
    library(package = "sfsmisc")
  4. Im Anschluss plotten wir die Daten erneut, aber ohne dabei die Achsen zu beschriften (xaxt="n", yaxt="n"). Die Beschriftung der Achsen wird im Anschluss mit Hilfe des Befehls eaxis() aus dem Paket SFSmisc hinzugefügt:
    plot(x, y, log="xy", xaxt="n", yaxt="n")
    sfsmisc::eaxis(side=1)   # X-Achse
    sfsmisc::eaxis(side=2)   # Y-Achse

    Im neuen Plot sind nun auf jedenfall acht Striche zur Unterteilung einer Größenordnung zu sehen:

    Erstellt man die Achsen eines Plots mit Hilfe des Befehls eaxis() aus dem Paket SFSmisc, so wird der Bereich zwischen zwei Größenordnungen in neun Teile unterteilt. Eventuell werden aber nicht nur die Zahlen für die Größenordnungen eingetragen, sondern auch für dazwischen liegende Werte (hier z.B. 5, 50 und 500).
    Erstellt man die Achsen eines Plots mit Hilfe des Befehls eaxis() aus dem Paket SFSmisc, so wird der Bereich zwischen zwei Größenordnungen in neun Teile unterteilt. Eventuell werden aber nicht nur die Zahlen für die Größenordnungen eingetragen, sondern auch für dazwischen liegende Werte (hier z.B. 5, 50 und 500).

    In diesem Beispiel wurde nicht explizit angegeben, welche Zahlen auf den Achsen zu sehen sein sollen. Daher sind nicht nur die Größenordnungen, sondern auch Werte dazwischen, als Zahlen eingetragen.

  5. Daher kann man nun noch explizit angeben, welche Zahlen sichtbar sein sollen (at=)
    plot(x, y, log="xy", xaxt="n", yaxt="n")
    sfsmisc::eaxis(side=1, at=10^c(0:3))
    sfsmisc::eaxis(side=2, at=10^c(-4,-2,0))

    In diesem Plot wurden nur Zahlen für die Größenordnungen angegeben. Auf der X-Achse wurde jede, auf der Y-Achse nur jede zweite Größenordnung mit einer Zahl versehen.
    In diesem Plot wurden nur Zahlen für die Größenordnungen angegeben. Auf der X-Achse wurde jede, auf der Y-Achse nur jede zweite Größenordnung mit einer Zahl versehen.

Transluzente Polygone in R-Plots einzeichnen

Manchmal möchte man Bereiche einer Grafik (z.B. einen Meßbereich) farblich hervorheben. Das kann man z.B. durch Polygone erreichen. Wichtig ist dabei, dass keiner der Datenpunkte vollständig überdeckt wird. Entweder zeichnet man die Polygone vor dem plotten der Datenpunkte ein, oder man verwendet teilweise durchsichtige Polygone (transluzente Polygone), durch die die Datenpunkte immer noch zu sehen sind.

nerv_test

Diese Transluzenz erreicht man in R dadurch, dass man einen achtstelligen hexadezimalen RGB-Farbwert verwendet. Die ersten 6 Stellen sind wie gewohnt für die drei Farben Rot, Grün und Blau reserviert. Die letzten beiden Stellen bestimmen die Deckfraft der zuvor angegebenen Farbe. Der Wert „00“ steht dabei für „vollständig durchsichtig (Farbe unsichtbar)“ und „ff“ steht für „100% deckend (undurchsichtig)“.

Im folgenden Beispiel werden Daten (x, y) geplottet.

x = 1:100
y = x + runif(n=length(x), min=-5, max=+5)
# Simulierte Standardabweichung:
sd = 10 + runif(n=length(x), min=-1, max=+1)
plot(x,y, las=1)

Die Standardabweichung (sd) soll in diesem Beispiel nicht als Fehlerbalken, sondern als farbiger Bereich um die Kurve geplottet werden. Das ist besonders dann hilfreich, wenn man eine hohe Dichte an Meßpunkten hat. Dazu werden der Funktion polygon() die x und y-Koordinaten übergeben. Für den Bereich „zuzüglich Standardabweichung“ (oberhalb der Messwerte; y+sd) gehen wir vom kleinsten zum größten x-Wert (x); für den Bereich „abzüglich Standardabweichung“ (unterhalb der Messwerte; y-sd) gehen wir vom größten zum kleinsten x-Wert (rv(x), und damit auch rev(y-sd)). Wird beschreiben das Polygon also in einer Rechtskurve.

# Bereiche der simulierten Standardabweichung einzeichnen
polygon(
  x = c(x, rev(x))
  , y = c(y+sd, rev(y-sd))
  , col = "#cc000033"
)

Auf diese Weise kann man auch einfache Formen, wie z.B. ein Rechteck, einfügen:

# Rechteck einzeichnen
polygon(
  x = c(40, 40, 60, 60)
  , y = c(-10, 200, 200, -10)
  , col = "#0000cc33"
)

Regressionsgerade in einen doppelt-logarithmischen Plot eintragen

Wer Potenzgesetze (PowerLaws) in Daten entdecken möchte, bestimmt die Steigung in einem Graphen mit logarithmischen Achsen (log-log-Plot). In der Regel möchte man diese Steigung dann in der Graphik als Regressionsgerade darstellen. In R ist genau diese Darstellung nicht ganz einfach, da die Funktion abline() hier versagt. Im folgenden Artikel möchte ich zeigen, wie man die Schätzung seines Modells dennoch in den Plot eintragen kann.

Wir haben Daten, die einem Potenzgesetz folgen (mit etwas statistischer Schwankung).

x = c(1:150)
y = x^-.5 * 155 + (runif(length(x), min=-3, max=3))

Diese Daten plotten wir in einer Graphik mit zwei logarithmischen Achsen.

plot(x, y, log="xy", cex=.5)

Da wir ein Potenzgesetz vermuten, berechnen wir die Parameter (Steigung und Y-Achsenabschnitt) des linearen Modells (lm()) der logarithmierten Daten (siehe auch Stackoverflow):

model = lm(log(y) ~ log(x))
model

In einem doppelt-logarithmischen Plot gibt es natürlich keinen Y-Achsenabschnitt (die Null wird nie erreicht). Beim Aufruf der Funktion abline(model) wird die Regressionsgerade (besonders bei verschobenen Funktionen) an der falschen Stelle dargestellt. Daher muss man für sein lineares Modell Werte vorhersagen (mit predict() bzw. predict.lm), die man im Anschluss in den Exponenten nimmt und als Linie (lines()) zum Plot hinzufügt (siehe auch Stackoverflow: Beitrag 1 und Beitrag 2).

# Schätzung für 2 Punkte machen
neuX=c(1e-10,1e10)
lines(neuX, exp(predict(model, newdata=list(x=neuX))) ,col="blue", type="o", pch=2)
 
# Schätzung für Datenpunkte machen
lines(x, exp(predict(model, newdata=list(x=x))), col="red", type="o", pch=4)

nerv_test

Hinweise: Der Variablenname in der Liste von newdata muss genau der selbe sein, wie der Name der Variablen, die beim Aufruf des linearen Modells verwendet wurde.

  1. Hat man seine Daten z.B. in einem data.frame gespeichert, sollte man das lineare Modell mittels with() aufrufen, anstatt die Bereiche mit Hilfe des $-Selektors auszuwählen.
  2. Wenn man das lineare Modell nur für einen Teil seiner Daten verwenden möchte, so sollte man die Option subset= von lm() verwenden, anstatt Teile des Vectors mit Hilfe der eckigen Klammern auszuwählen.

Die folgenden Beispiele zeigen zwei Mögliche Schreibweisen, eine lineare Regression für die die Datensätze 40 bis 80 einzuzeichnen. Die Daten dafür sind:

1
2
3
4
x = c(1:150)
y = x^-.5 * 155 + (runif(length(x), min=-3, max=3))
daten = data.frame(x,y)
range = c(40:80)

Beispiel 1: with() wird nur für den Aufruf des lm() verwendet. die Ergebnisse der Regression werden in model gespeichert.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
plot(daten$x, daten$y, log="xy", cex=.5)
model = with(
  data = daten,
  expr = { model = lm(log(y) ~ log(x), subset=range); model }
)
model
 
lines(
  x[range]
  , exp(predict(model, newdata=list(x=x[range])))
  , col="red"
  , type="o"
  , pch=4
)

Beispiel 2: with() wird für die gesamte Behandlung der Daten verwendet (plot(), lm(), lines()).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
with(
  data = daten,
  expr = {
    plot(x, y, log="xy", cex=.5)
    model = lm(log(y) ~ log(x), subset=range)
    lines(
      x[range]
      , exp(predict(model, newdata=list(x=x[range])))
      , col="red"
      , type="o"
      , pch=4
    )
  }
)

Eine Liste zu einem data.frame konvertieren in R

Manchmal liegen in R Daten in Form einer Liste vor (z.B. die Ausgabe der Funktion hist()). In einer Liste können Daten enthalten sein, die eine unterschiedliche Länge haben. Im Falle von hist() sind z.B. Arrays unterschiedlicher Länge gemeinsam mit einfachen Variablen in einer Liste kombiniert.

Für weitere Analysen kann es jedeoch hilfreich sein, die gleich langen Teile in einen separaten data.frame auszulagern. Dann kann z.B. mit Hilfe der Funktion subset() relativ einfach Bereiche aus den Daten herausfiltern.

In R gibt es keine vorgefertigte Funktion, die direkt Daten vom Typ list in Daten vom Typ data.frame konvertieren kann. Das liegt daran, dass in einem data.frame alle enthaltenen Variablen (in der Regel Arrays) die gleiche Länge haben müssen. Unsere Entscheidung, welche Variablen aus der Liste in den data.frame übernommen werden sollen, müssen wir R mitteilen.

Das folgende Beispiel zeigt, wie man das sehr einfach bewerkstelligen kann:

1
2
3
4
5
6
7
8
# Liste mit Arrays ("x", "y", "sd") und einer Variablen ("name")
liste = list(x=c(1:5), y=c(11:15), sd=c(0,1,2,0,0), name="Meine Liste")
 
# data.frame ("df") erzeugen mit den Variablen "x", "y" und "sd"
df = with(data=liste, expr=data.frame(x, y))
 
# Struktur des data.frames ausgeben
str(df)
'data.frame':	5 obs. of  3 variables:
 $ x : int  1 2 3 4 5
 $ y : int  11 12 13 14 15
 $ sd: num  0 1 2 0 0

Im Anschluss kann man sehr einfach nur bestimmte Teile des data.frame ausgeben lassen.
Das kann man z.B. direkt über den Index (df[zeile, spalte]) oder über die Funktion subset() machen.

1
2
3
4
5
6
7
8
9
# Nur das zweite und vierte Element (Zeile) anzeigen
df[c(2,4),]
 
# Nur die Elemente (Zeilen) anzeigen, bei denen x>2 ist
df[(df$x>2),]
 
# Nur die Elemente (Zeilen) der Variablen (Spalten) x und y ausgeben,
# bei denen sd<1 und x>2 ist
subset(df, subset=(x>2 & sd<1), select=c(x,y))

Klassifizierung von Messwerten in R

Stell Dir vor, Du hast die Länge von 1000 Fischen gemessen. Im Anschluss möchtest Du die eine Häufigkeitsverteilung (Histogramm) der Größen erstellen. Je nachdem wie genau du gemessen hast, wirst du keine zwei Fische mit der gleichen Länge finden. Daher bist Du gut beraten, die Daten zunächst in bestimmte Längenklassen einzuteilen (z.B. „Anzahl von Fischen zwischen 23cm und 24cm“). Für diese Klassifizierung (binning) steht Dir in R die Funktion hist() zur Verfügung.

Nehmen wir mal an, die Längen der Fische folgen einer Normalverteilung. Im Durschnitt haben die Fische eine Länge von 25cm (± 5cm)

1
2
3
# Ziehe Eintausend Zufallszahlen aus einer Normalverteilung
# (Mittelwert: 25; Standardabweichung: 5)
laengen = rnorm(n=1e3, mean=25, sd=5)

Mit der Funktion hist() kannst Du die Daten nun in Klassen einteilen und plotten lassen.

1
2
3
# Klassifiziere die Daten
# (=Erstelle eine Histogramm und stelle es dar)
gebinnt = hist(laengen, plot=TRUE)
Automatisch erstelltes Histogramm der Beispieldaten. Die Klassengrenzen wurden von R bestimmt.
Automatisch erstelltes Histogramm der Beispieldaten. Die Klassengrenzen wurden von R bestimmt.

hist() erstellt nun eine list, in der die Klassengrenzen (breaks), die Häufigkeiten (counts), Dichten (densitiy) und Klassenmitten (mids), sowie der Name der ursprünglichen Variable (xname) und die Information, ob die Klassen alle gleich groß sind (equidist), gespeichert werden:

1
str(gebinnt)
List of 7
 $ breaks     : num [1:9] 5 10 15 20 25 30 35 40 45
 $ counts     : int [1:8] 3 29 140 370 312 129 15 2
 $ intensities: num [1:8] 0.0006 0.0058 0.028 0.074 0.0624 0.0258 0.003 0.0004
 $ density    : num [1:8] 0.0006 0.0058 0.028 0.074 0.0624 0.0258 0.003 0.0004
 $ mids       : num [1:8] 7.5 12.5 17.5 22.5 27.5 32.5 37.5 42.5
 $ xname      : chr "laengen"
 $ equidist   : logi TRUE
 - attr(*, "class")= chr "histogram"

Für weitere Analysen kannst Du selbst noch die Breite der Klassen (breite) und den Anteil der Klassen an der Stichprobe (anteil) zur Liste hinzufügen:

1
2
gebinnt$breite = diff(gebinnt$breaks)
gebinnt$anteil = gebinnt$counts / sum(gebinnt$counts)

Eigene Klassenbreiten festlegen

Du hast natürlich auch die Möglichkeit, selbst Klassengrenzen (breaks) anzugeben. Die Klassen müssen auch nicht alle gleich Groß sein. Bei einigen Daten können z.B. logarithmische Bin-Größen sinnvoll sein.

1
2
3
4
# Klassen mit Klassenbreite 1 (cm) erstellen
gebinnt = hist(
  laengen, plot=T, breaks=c(0:50)
  , xlab="Länge (cm)", ylab="Häufigkeit", main="Häufigkeitsverteilung")
Histogramm mit selbst gewählten Klassengrenzen (hier: Klassenbreite=1cm)
Histogramm mit selbst gewählten Klassengrenzen (hier: Klassenbreite=1cm)

Hierbei muss man allerdings darauf achten, dass keiner der Werte ausgeschlossen wurde, weil er ausserhalb der gewählten Klassen lag bzw. auf die untere (oder obere) Klassengrenze gefallen ist.

Fehler in hist.default(laengen, plot = T, breaks = c(15:50), xlab = "Länge (cm)",  : 
  einige 'x' nicht gezählt: evtl. überdecken die 'breaks' nicht den gesamten Bereich von 'x'

Daher sollte man immer prüfen, ob die Summe der Werte in den Klassen auch tatsächlich der Stichprobengröße entspricht.

# Sind alle Messwerte der Stichpobe im Histogramm berücksichtigt?
sum(gebinnt$counts)

Eventuell muss man dann weitere Klasse hinzunehmen, bzw. die untere (oder obere) Klassengrenze zu einer Klasse hinzu zählen (Option include.lowest=T).

Alternativen

Alternativ zur Funktion hist() könnte man ein Histogramm auch durch die Anwendung der Funktion table() auf gerundete Daten (round()) erstellen: table(round(daten)).

Die Funktion table() ist dazu gedacht die „Anzahl gleicher Werte zu ermitteln“. Das entspricht aber nicht ganz dem Gedanken der Klassifizierung.

Fehlerbalken in R-Grafik einzeichnen

In wissenschaftlichen Publikationen werden häufig Messwerte in Grafiken dargestellt. Messungen sind allerdings immer fehlerbehaftet. Daher sollte man dem Leser mitteilen, wie verlässlich die gemessenen Werte sind. In der deskriptiven Statistik gibt es dafür diverse Kennzahlen, wie z.B. die Standardabweichung oder den Standardfehler. Diese werden oft als Fehlerbalken in die Grafiken eingefügt.

errorbars
Grafisch dargestellte Messwerte mit Fehlerbalken.

In R gibt es leider keine Standardmethode, die diese Aufgabe übernimmt. Zunächst hatte ich nach Anleitungen im Internet eigene Funktionen für diese Aufgabe erstellt.

  • Eine Lösung fand ich bei Stackoverflow. Allerdings wurden hier die horizontalen Linien bei einer logarithmischen X-Achse nicht korrekt angezeigt (links und rechts waren die Linien unterschiedlich lang).
  • Die Lösung von MonkeysUncle war in der Hinsicht besser. Allerdings wurde keine horizontale Linie angezeigt, wenn der Fehler 0 war. Zudem wurde gar kein Fehlerbalken angezeigt, wenn die untere Grenze des Fehlers (bei logarithmischer Y-Achse) im negativen Bereich lag.

Schließlich habe ich die Funktion errbar() im Paket Hmisc gefunden (siehe auch StackOverflow). Hier werden:

  1. Die horizontalen Linien auch bei logarithmischer X-Achse korrekt gezeichnet.
  2. Bei logarithmischer Y-Achse und negativer unteren Fehlergrenze, wird zumindest die obere Fehlergrenze geplottet.
  3. Eine horizontale Linie wird auch hinzugefügt, wenn der Fehler 0 ist.
  4. Die Y-Ranges werden automatisch gewählt, wenn die Fehlerbalken ausserhalb der Grafik liegen würden.

Das folgende Beispiel zeigt, wie man Datenpunkte mit ihren Fehlerbalken in eine Grafik einträgt.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Datenpunkte (x, y) mit einer Standardabweichung (sd)
daten = data.frame(
  x  = c(1:5)
  , y  = c(1.1, 1.5, 2.9, 3.8, 5.2)
  , sd = c(0.2, 0.3, 0.2, 0.0, 0.4)
)
 
# Paket "Hmisc" installieren und laden 
install.packages("Hmisc", dependencies=T)
library("Hmisc")
 
# Plot mit automatischer Wahl des Y ranges
with (
  data = daten
  , expr = Hmisc::errbar(x, y, y+sd, y-sd, pch=1)
)
 
# Datenpunkte mit Fehlerbalken zu einem
# existierenden Koordinatensystem hinzufügen (add=T)
plot(daten$x, daten$y, type="n", xlab="X-Werte", ylab="Y-Werte")
with (
  data = daten
  , expr = Hmisc::errbar(x, y, y+sd, y-sd, pch=2, add=T)
)

Die Funktion errbar() erweitert die Funktion plot(). Daher können die Grafikparameter von plot() und par() (hier z.B. pch) verwendet werden. Zusätzlich gibt es noch Parameter, mit denen sich z.B. die breite der horizontalen Begrenzung (cap) steuern lässt.

Funktion aus einem bestimmten R-Paket laden

Die Programmiersprache R kann durch Funktionen erweitert werden, die von anderen Nutzern erstellt wurden. Oftmals sind diese Erweiterungen in Paketen zusammengefasst. Ein Problem tritt allerdings dann auf, wenn in unterschiedlichen Paketen Funktionen enthalten sind, die den gleichen Namen tragen (siehe: Stackoverflow).

Ich nutze z.B. folgende Pakete sehr häufig:

  • sfsmisc – Mit der Funktion eaxis() lassen sich schöne logarithmische Achsen erstellen. Die Beschriftung wird dabei von der Funktion pretty10exp() übernommen.
  • Hmisc– Mit der Funktion errbar() lassen sich sehr leicht Fehlerbalken in Plots eintragen.
  • TeachingDemos – Mit der Funktion subplot() lassen sich sehr einfach Insets zu Plots hinzufügen.

Lädt man allerdings die Pakete,

## Installation der Pakete
install.packages("sfsmisc", dependencies=T)
install.packages("Hmisc", dependencies=T)
install.packages("TeachingDemos", dependencies=T)
 
## Laden der Pakete
library(package = "sfsmisc")
library(package = "Hmisc")
library(package = "TeachingDemos")

dann erhält man folgende Fehlermeldungen:

Attache Paket: ‘Hmisc’

Das folgende Objekt ist maskiert from ‘package:survival’:
    untangle.specials

Das folgende Objekt ist maskiert from ‘package:sfsmisc’:
    errbar

Das folgende Objekt ist maskiert from ‘package:base’:
    format.pval, round.POSIXt, trunc.POSIXt, units

Man kann erkennen, dass sowohl Hmisc als auch sfsmisc eine Funktion namens errbar() enthalten.

Attache Paket: ‘TeachingDemos’

Das folgende Objekt ist maskiert from ‘package:Hmisc’:
    cnvrt.coords, subplot

Man kann erkennen, dass sowohl TeachingDemos als auch Hmisc eine Funktion mit Namen subplot() enthalten.

Um nun die Funktion aus einem bestimmten Paket zu verwenden, kann man R dies mit Hilfe von zwei Doppelpunkten :: mitteilen:

# Verwende die Funktion subplot() aus dem Paket TeachingDemos
TeachingDemos::subplot(...)
 
# Verwende die Funktion errbar() aus dem Paket Hmisc
Hmisc::errbar(...)

Meine Empfehlung: Wenn man Funktionen aus einem Paket verwendet, sollte man beim Aufruf der Funktion immer diese Notation verwenden. Zum einen hilft es dem Leser zu erkennen, woher diese Funktion stammt. Zum anderen verhindert man unnötige Verwirrung, wenn man zu einem späteren Zeitpunkt Pakete einbindet, die Funktionen mit gleichem Namen enthalten.