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)

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.

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
    )
  }
)