Wie wir gesehen haben, funktioniert unsere selbst erstellte Funktion bereits einwandfrei. Nach dem Prinzip der Sparsamkeit – und zum Teil auch aus stilistischen Gründen – lässt sich die Schreibweise jedoch noch weiter optimieren.

Ein guter Ausgangspunkt dafür ist die Bedingung in der runden Klammer der if-Abfrage. Wir haben gelernt: Wird die Bedingung als TRUE ausgewertet, wird der Code in den geschweiften Klammern ausgeführt – wir die Bedingung als FALSE ausgewertet, wird zum nächsten Abschnitt, etwa einem else, weitergesprungen. In unserem bisherigen Beispiel prüfen wir mit einer logischen Abfrage, ob empirical == TRUE ist. Diese Abfrage liefert entweder TRUE oder FALSE zurück.

Doch warum ist das überhaupt nötig? Das Objekt empirical steht ja bereits selbst für einen logischen Wert – entweder TRUE oder FALSE. Eine zusätzliche Abfrage ist also überflüssig. Statt if (empirical == TRUE) reicht schlicht if (empirical).

var_eigen <- function(x, empirical){
  n <- length(x)
  x_quer <- mean(x)
  if (empirical) {
    var <- sum((x - x_quer)^2) / n
  } else {
    var <- sum((x - x_quer)^2) / (n - 1)
  }
  return(var)
}

Außerdem fällt auf, dass unter beiden Bedingungen – also im if- und im else-Zweig – ein Großteil der Berechnung identisch ist. Bereits zu Beginn dieses Tutorials haben wir gelernt, dass sich die empirische Varianz und die Populationsvarianz sehr einfach ineinander überführen lassen. Daher stellt sich die Frage: Brauchen wir überhaupt eine vollständige Aufteilung in if und else? Eleganter wäre es, zunächst grundsätzlich die empirische Varianz zu berechnen – unabhängig von der Einstellung des Arguments empirical. Nur falls empirical == FALSE, wenden wir anschließend die Korrektur für die Populationsvarianz an, indem wir die Berechnung entsprechend anpassen. Das macht den Code nicht nur kürzer, sondern auch deutlich lesbarer.

var_eigen <- function(x, empirical){
  n <- length(x)
  x_quer <- mean(x)
  var <- sum((x - x_quer)^2) / n
  if (!empirical) {
    var <- var * (n / (n - 1))
  }
  return(var)
}

Dabei fällt auf, dass der Code in der if-Bedingung nur dann ausgeführt werden soll, wenn empirical auf FALSE gesetzt ist. Wie wir bereits gelernt haben, lässt sich ein logischer Wert in sein Gegenteil umkehren, indem man ein Ausrufezeichen (!) davor setzt. So wird aus FALSE ein TRUE – und genau das benötigen wir hier: Mit if (!empirical) stellen wir sicher, dass der Korrektur-Schritt nur dann ausgeführt wird, wenn keine empirische, sondern die Populationsvarianz berechnet werden soll.

Im weiteren Verlauf des Tutorials arbeiten wir jedoch wieder mit der ursprünglich erstellten Version der Funktion, die ohne in dieser Vertiefung vorgestellten Optimierungen auskommt. Zur Erinnerung, hier noch einmal der vorherige Stand:

var_eigen <- function(x, empirical){
  n <- length(x)
  x_quer <- mean(x)
  if (empirical == TRUE) {
    var <- sum((x - x_quer)^2) / n
  } else {
    var <- sum((x - x_quer)^2) / (n - 1)
  }
  return(var)
}

Hinweis: Der folgende Abschnitt ist unabhängig von der vorherigen Vertiefung zur Optimierung der Funktion lesbar und verständlich, kann aber selbstverständlich mit dieser kombiniert werden.

Das Vergessen der Eingabe eines Arguments ist nur eine von vielen möglichen Fehlerquellen bei der Verwendung unserer Funktion – und nicht immer werden solche Fehler sofort deutlich. So könnte es beispielsweise passieren, dass Nutzer:innen dem Argument empirical einen Wert zuweisen, den wir gar nicht vorgesehen haben (etwa einen Text oder eine Zahl oder einen Text statt TRUE oder FALSE). Testen wir das:

var_eigen(x = fb24$time_pre, empirical = "ja")

## [1] 425.3645

var_eigen(x = fb24$time_pre, empirical = 200)

## [1] 425.3645

Interessanterweise erhalten wir in solchen Fällen dennoch ein Ergebnis – was zunächst verwundern mag. Ein Blick in unseren Funktionscode zeigt jedoch, woran das liegt: In der if-Abfrage prüfen wir lediglich, ob empirical == TRUE ist. Das ist in den getesteten Fällen eindeutig nicht gegeben, weshalb automatisch der else-Zweig ausgeführt wird. Unsere Funktion läuft also durch, erzielt aber nicht den beabsichtigten Effekt.

Um solche unerwünschten Effekte zu vermeiden, können wir eine eigene Fehlermeldung definieren, die ausgegeben wird, wenn ein Argument nicht den erwarteten Wert hat. Dazu verwenden wir die Funktion stop(). Diese wird typischerweise in einer if-Abfrage eingesetzt, um beispielsweise zu prüfen, ob ein Argument vom richtigen Typ ist. In unserem Fall können wir testen, ob der Eintrag in empirical tatsächlich ein logischer Wert (logical) ist.

var_eigen <- function(x, empirical = TRUE){
  if(is.logical(empirical) == FALSE){
    stop("Das Argument 'empirical' muss TRUE oder FALSE sein.")
    }
  n <- length(x)
  x_quer <- mean(x)
  if (empirical == TRUE) {
    var <- sum((x - x_quer)^2) / n
  } else {
    var <- sum((x - x_quer)^2) / (n - 1)
  }
  return(var)
}

Im weiteren Verlauf des Tutorials arbeiten wir jedoch wieder mit der ursprünglich erstellten Version der Funktion, die ohne in dieser Vertiefung vorgestellten Optimierungen auskommt. Zur Erinnerung, hier noch einmal der vorherige Stand:

var_eigen <- function(x, empirical = TRUE){
  n <- length(x)
  x_quer <- mean(x)
  if (empirical == TRUE) {
    var <- sum((x - x_quer)^2) / n
  } else {
    var <- sum((x - x_quer)^2) / (n - 1)
  }
  return(var)
}

Logische Abfragen und Bedingungen: if und else

Im Tutorial haben wir bereits gesehen, dass sich mit if und else sogenannte Wenn-Dann-Bedingungen in den Code einbauen lassen. An dieser Stelle wollen wir uns dem Thema noch einmal aus einer etwas technischeren Perspektive widmen. Die folgenden Beispiele mögen teilweise etwas konstruiert wirken – wir haben jedoch versucht, sie möglichst anschaulich zu gestalten. Für alle, die sich intensiver mit R beschäftigen möchten – etwa, um eigene Funktionen zu schreiben oder bestehende Funktionen näher zu inspizieren – ist das folgende Wissen auf jeden Fall hilfreich und essenziell.

if-Abfragen

Einfache if-Abfrage

Im Tutorial haben wir – abgesehen von den Vertiefungen – if und else nur als gemeinsame Befehle verwendet. Dabei kann if auch eigenständig eingesetzt werden, denn es bedeutet lediglich, dass eine bestimmte Konsequenz nur dann ausgeführt wird, wenn die zugehörige Bedingung erfüllt ist.

a <- 3 # Zunächst definieren wir eine Variable
# Nutzung einer einfachen if-Abfrage
if (a == 3) {
  print("Ja, die Variable a enthält den Wert 3")
}

## [1] "Ja, die Variable a enthält den Wert 3"

Für das Verständnis solcher Abfragen ist es hilfreich, die verschiedenen Schritte der Syntax einzeln zu betrachten. Das gilt auch für die restlichen Themen dieses Blocks. In R-Studio können Sie einzelne Abschnitte des Codes markieren und ausführen, um zu testen, was diese beinhalten. Füren Sie Ihren Code immer in kleinen Schnipseln aus, wenn Sie verstehen wollen, was passiert. Hier wird im ersten Schritt die Bedingung evaluiert:

(a == 3)

## [1] TRUE

In diesem Fall stimmt die logische Abfrage (a enthält tatsächlich den Wert 3), und wir erhalten in der Konsole den Output TRUE. Deshalb wird der danach definierte print-Befehl ausgeführt.

Wenn jedoch a einen anderen Wert enthält, trifft die Bedingung nicht zu (FALSE) und der folgende Befehl wird deshalb nicht ausgeführt.

a <- 5
if (a == 3) {
  print("Ja, die Variable a enthält den Wert 3")
}

Da die Bedingung nicht erfüllt ist, wird die Konsequenz nicht ausgeführt – beim Ausführen erscheint daher lediglich der Befehl in der Konsole, jedoch kein Ergebnis.

if-Abfrage mit mehreren Möglichkeiten

Oft haben wir aber mehrere Argumente, die untersucht werden können. Bspw. können wir testen, ob ein Wert sich in einer Liste wiederfindet. Wenn wir beispielsweise herausfinden wollen, ob die Person, die in der Variable person gespeichert ist, ein Hauptcharakter aus der Serie Friends ist, können wir dies mit dem folgenden Befehl tun:

person <- "Monica"
if (person %in%  c("Monica", "Rachel", "Chandler",  "Phoebe", "Ross", "Joey")) {
  print("Yes, this is a character from Friends.")
}

## [1] "Yes, this is a character from Friends."

Hier erhalten wir die Antwort, ja, Monica ist eine Figur aus der Serie. Der Ausdruck %in% steht sinngemäß für “ist ein Element aus der folgenden Auswahl”.

Wenn wir die gleiche Abfrage auf eine andere Person anwenden, trifft die Bedingung nicht zu, und der Befehl wird nicht ausgeführt (in der Konsole erscheint also nichts außer der Befehl).

person = c("Marcus")
if (person %in%  c("Monica", "Rachel", "Chandler",  "Phoebe", "Ross", "Joey")) {
  print("Yes, this is a character from Friends.")
}

Genauso könnten wir aber auch eine Liste von Personen haben und uns entweder fragen, ob mindestens eine Person aus dieser Liste bei “Friends” mitgewirkt hat, oder ob alle Personen dort mitgewirkt haben. Dies geht mit any oder all:

persons <- c("Monica", "Marcus")
if (any(persons %in%  c("Monica", "Rachel", "Chandler",  "Phoebe", "Ross", "Joey"))) {
  print("Yes, at least one of them is a character from Friends.")
}

## [1] "Yes, at least one of them is a character from Friends."

if (all(persons %in%  c("Monica", "Rachel", "Chandler",  "Phoebe", "Ross", "Joey"))) {
  print("Yes, all of them are a character from Friends.")
}

Schauen wir uns nur die logische Abfrage an, um den Unterschied der beiden Befehle deutlich zu machen.

any(persons %in%  c("Monica", "Rachel", "Chandler",  "Phoebe", "Ross", "Joey")) 

## [1] TRUE

# mindestens 1 ist TRUE

all(persons %in%  c("Monica", "Rachel", "Chandler",  "Phoebe", "Ross", "Joey"))

## [1] FALSE

# alle sind TRUE

any liefert in einer logischen Abfrage den Wert TRUE, wenn mindestens ein Element die Bedingung erfüllt, während all nur dann TRUE ergibt, wenn alle Elemente die Bedingung erfüllen. Mit diesen beiden Befehlen lässt sich also leicht prüfen, ob mindestens ein oder alle Elemente eine bestimmte Eigenschaft aufweisen.

Abgleich mit einem Datum

Es sind beispielsweise auch logische Abfragen mit Zeitpunkten und Daten möglich. Zum Beispiel können wir mit dem Befehl weekdays(Sys.Date()) ermitteln, welcher Wochentag gerade ist, und dann abgleichen, ob Freitag ist. Wenn Sie diesen Befehl selbst testen, achten Sie darauf, ob nach der Voreinstellung Ihres Rechners der aktuelle Wochentag auf Englisch oder Deutsch ausgegeben wird. Das finden Sie heraus, indem Sie nur den kleinen Codeabschnitt weekdays(Sys.Date()) ausführen.

if (weekdays(Sys.Date()) == "Friday") {
  print("Fast Wochenende!")
}

Verknüpfung logischer Abfragen

Wie im letzten Semester bereits besprochen, können logische Bedingungen mit & (logisches “und”) und | (logisches “oder”) verknüpft werden. Wenn die gesamte logische Abfrage als Ergebnis TRUE zurückgibt, wird die R-Syntax in den geschwungenen Klammern ausgeführt; wenn es FALSE ergibt, passiert nichts. Zum Beispiel könnten wir so testen, ob entweder Samstag oder Sonntag ist und herausfinden, ob wir uns freuen dürfen.

if (weekdays(Sys.Date()) == "Saturday" | weekdays(Sys.Date()) == "Sunday") {
  print("Hoch die Hände, Wochenende!")
}

Durch die logische Verknüpfung mit | (logisches “oder”) wird die gesamte Abfrage wahr, wenn entweder der erste oder der zweite Teil zutrifft (in Worten “ist heute entweder Samstag oder Sonntag?”). Hätten wir an dieser Stelle stattdessen eine Verknüpfung mit & (logisches “und”) gewählt, könnte die Bedingung “ist heute Samstag und Sonntag?” nicht zutreffen, würde also immer FALSE zurückgeben. An anderer Stelle ist das & aber notwendig, wenn mehrere Bedingungen erfüllt sein sollen.

Bei der Verknüpfung dieser logischen Abfragen muss auf Klammersetzung geachtet werden, wenn die Verknüpfung komplizierter wird. Beispiel: “Ist heute (Samstag oder Sonntag) und scheint die Sonne?”. Als Übung können Sie versuchen diese logische Abfrage mit fiktiven Variablen in Code auszudrücken.

Abgleich mit mehreren Alternativen: if-else-Abfragen

Häufig wollen wir nicht nur konditional einen Befehl ausführen, oder nicht ausführen, sondern möchten einen anderen Befehl angeben, der ausgeführt wird, wenn die Bedingung nicht zutrifft (wie im Tutorial für die Erstellung der Funktion verwendet). Wir ergänzen also das else, das zum Tragen kommt, wenn die Bedingung nicht zutrifft. Dies lässt sich fast wörtlich lesen “If the condition is true, then do one thing. Otherwise (else), do the other thing.”

# mehrere Zeilen
if (weekdays(Sys.Date()) == "Saturday" | weekdays(Sys.Date()) == "Sunday") {
  print("Hoch die Hände, Wochenende!")
} else {
  print("Nur noch wenige Tage bis zum ersehnten Wochenende!")
}

## [1] "Nur noch wenige Tage bis zum ersehnten Wochenende!"

Bei Code über mehrere Zeilen ist es wichtig, die geschweiften Klammern korrekt zu setzen. Nach der Bedingungsabfrage öffnen sich geschweifte Klammern, die den ersten konditionalen Befehl einschließen. Das else folgt darauf. Danach wird der alternative Befehl wieder in geschweiften Klammern eingefasst. Der else-Befehl muss, wenn wir if und else in reinem Skript verwenden, in der gleichen Zeile stehen wie die geschlossene geschweifte Klammer (und auch in einer Funktion ist es empfehlenswert).

else if-Bedingungen

Häufig werden mehrere Abfragen ineinander geschachtelt, sodass die Ausdrücke schnell sehr kompliziert werden können. Falls in mehreren Schritten verschiedene Bedingungen abgefragt werden, und verschiedene Konsequenzen folgen sollen, kann auch das else if verwendet werden. Hierbei werden verschiedene Möglichkeiten abgefragt, für die verschiedene Befehle ausgeführt werden sollen. Wenn die erste Bedingung nicht zutrifft, wird die zweite Bedingung (nach dem else if) geprüft, wenn diese auch nicht zutrifft, wird das nächste else if geprüft. Der Befehl nach dem else wird dann nur ausgeführt, wenn keine der vorherigen Bedingungen zutrifft.

Hier sehen Sie ein Beispiel für eine if-else-Abfrage, die Sie jeden Morgen nutzen können, um herauszufinden, wie Sie sich heute fühlen sollten.

if (weekdays(Sys.Date()) %in% c('Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday')) {
  if (weekdays(Sys.time()) == 'Monday') {
    print('Zurück ins Bett...')
    } else if (weekdays(Sys.time()) == 'Wednesday') { 
    print('Wuhu, es ist Mitte der Woche!')
      } else if (weekdays(Sys.time()) == 'Friday') { 
    print('Yeah, das Wochenende steht bevor!')
        } else {
    print('Es ist irgendein anderer Tag.')
        }
  } else {
  print("Hoch die Hände, Wochenende!")
}

## [1] "Hoch die Hände, Wochenende!"

Wir versuchen nachzuvollziehen, was in dieser verschachtelten if-else-Abfrage passiert. Zunächst wird geprüft, ob es sich heute um einen Wochentag handelt:

if (weekdays(Sys.Date()) %in% c('Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday')) {

Wenn dem so ist, wird der nächste Block ausgeführt:

if (weekdays(Sys.time()) == 'Monday') {
    print('Zurück ins Bett...')
    } else if (weekdays(Sys.time()) == 'Wednesday') { 
    print('Wuhu, es ist Mitte der Woche!')
      } else if (weekdays(Sys.time()) == 'Friday') { 
    print('Yeah, das Wochenende steht bevor!')
        } else {
    print('Es ist irgendein anderer Tag.')
        }

Dieser fragt ab, ob heute Montag ist (if (weekdays(Sys.time()) == 'Monday') {), falls dem nicht so ist, wird geprüft ob Mittwoch ist (} else if (weekdays(Sys.time()) == 'Wednesday') {), wenn dem wiederum nicht so ist, dann wird geprüft ob Freitag ist (} else if (weekdays(Sys.time()) == 'Friday') {) und falls dem auch nicht so ist, dann wird else ausgeführt.

Ist heute kein Wochentag, dann wird direkt die else Bedingung für die erste if-Abfrage ausgedruckt, nämlich

  } else {
  print("Hoch die Hände, Wochenende!")
}

Dieses Spiel der Verschachtelung lässt sich beliebig erweitern.

Funktion ifelse

Wenn nur eine Bedingung abgefragt werden soll, und je nach Ergebnis einer von zwei Befehlen folgen soll, kann der Code abgekürzt werden. Für einzelne Ereignisse kann in R die Notation mithilfe der ifelse()-Funktion verwendet werden. Diese ist (anders als die if-Abfragen) eine klassische Funktion mit Argumenten. Die Funktion nimmt drei Argumente entgegen:

• test: die Bedingung
• yes: was getan werden soll, wenn die Bedingung zutrifft
• no: was getan werden soll, wenn die Bedingung nicht zutrifft

ifelse(test = weekdays(Sys.Date()) == 'Friday', yes = 'Yeah, das Wochenende steht bevor!', no = 'Es ist irgendein anderer Tag...')

## [1] "Es ist irgendein anderer Tag..."

So wird die gleiche if-else-Abfrage verkürzt dargestellt. Gerade für komplexere Abfragen kann aber nicht immer diese verkürzte Form gewählt werden. Die längere Version ist immer dann von Vorteil, wenn der auszuführende R-Code mehrere Zeilen lang ist oder z.B. weitere Bedingungen enthält. Auch lässt er sich oft übersichtlicher gestalten.

Schachtelung von for-Loops

for-Loops können auch ineinander geschachtelt werden. Dabei wird für die zweite Iteration häufig ii als Platzhalter verwendet. Im Befehl kann dann auf i und ii Bezug genommen werden. Hier sehen Sie beispielsweie, wie Sie ineinander geschachtelt durch einen Vektor aus Buchstaben und einen Vektor aus Zahlen iterieren. Was passiert, wenn Sie den ersten print-Befehl außerhalb des inneren Loops platzieren? Versuchen Sie, den Unterschied nachzuvollziehen.

Buchstaben <- c("A", "B", "C")
Zahlen <- c(1,2)
for (i in Buchstaben) {
  for (ii in Zahlen) {
    print(i)
    print(ii) 
  }
}

## [1] "A"
## [1] 1
## [1] "A"
## [1] 2
## [1] "B"
## [1] 1
## [1] "B"
## [1] 2
## [1] "C"
## [1] 1
## [1] "C"
## [1] 2

Das i und ii sind hier willkürlich gewählte Platzhalter. Wir könnten auch jeden anderen Buchstaben (oder Zeichenkombination) wählen. Es bietet sich an, Namen zu vergeben, die sinnvoll für den Sachverhalt sind und es auch bei mehreren Schachtelungen möglich machen, zu erkennen, um was es sich gerade handelt. Für dieses Beispiel würde sich bspw. for(buchstabe in Buchstaben) [oder auch: for(b in Buchstaben)] und for(zahl in Zahlen) [oder auch: for(z in Zahlen)] anbieten.

while-Loops

In while-Loops wird der Code so lange ausgeführt, bis eine vorab definierte Bedingung erfüllt ist. Ein einfaches Beispiel wäre, so lange einen Münzwurf zu simulieren, bis man 10 mal “Kopf” geworfen hat. Dafür müssen wir zum Einen die Münze als Objekt mit zwei Auswahlmöglichkeiten Kopf und Zahl anlegen, und ein leeres Objekt, in das wir die Ergebnisse der Münzwürfe speichern können.

# Münze erstellen
coin <- c('Kopf', 'Zahl')

# Leeres Objekt für die Aufzeichnung erstellen
toss <- NULL

Als nächstes schreiben wir den eigentlichen Loop. Dieser enthält eine logische Abfrage, die abfragt, ob die Anzahl der Kopf-Würfe unter 10 ist. Führen Sie nacheinander die Codeabschnitte toss == 'Kopf', sum(toss == 'Kopf') und sum(toss == 'Kopf')<10 aus, um zu verstehen, wie sich die logische Abfrage zusammensetzt. (Hinweis: den logischen Werten TRUE und FALSE sind die Zahlen 1 und 0 zugeordet.)

# Loop
while (sum(toss == 'Kopf')<10) {
  toss <- c(toss, sample(coin, 1))
}

# Würfe ansehen
toss

##  [1] "Zahl" "Zahl" "Kopf" "Zahl" "Zahl" "Zahl" "Kopf" "Kopf" "Zahl" "Kopf" "Zahl" "Zahl" "Kopf" "Zahl"
## [15] "Zahl" "Kopf" "Kopf" "Zahl" "Kopf" "Kopf" "Kopf"

repeat-Loops

Im Gegensatz zu for und while wird bei repeat zunächst kein explizites Abbruchkriterium definiert. Stattdessen wird repeat häufig genutzt, wenn es verschiedene oder veränderliche Abbruchkriterien für den Loop gibt. Diese Kriterien werden bei repeat allerdings innerhalb des Loops definiert - in den meisten Fällen wird dazu über if mindestens eine Bedingung definiert, unter der die Ausführung abgebrochen werden soll.

Ein einfaches Beispiel hierfür ist es, eine Fibonacci-Sequenz zu bilden (eine Sequenz in der eine Zahl immer die Summe der vorherigen beiden Zahlen ist):

$$a_n := a_{n-1} + a_{n-2}$$ für $n>1$ und $a_1=a_0=1$.

Wir können nun repeat nutzen, um die Sequenz abzubrechen, wenn die letzte Zahl z.B. größer als 1000 ist. An dieser Stelle wissen wir nicht, welches Element das sein wird, bzw. nach wie vielen Schritten dies passiert, wodurch es geschickter ist, innerhalb des Loops das Kriterium zu evaluieren. Wir nutzen hier n-1, n, und n+1 als Schritte, da es das 0-te Element in Vektoren in R nicht gibt.

fibo <- c(1, 1)

repeat {
  n <- length(fibo)
  fibo[n+1] <- fibo[n] + fibo[n - 1]
  if (fibo[n+1] > 1000) break
}

fibo

##  [1]    1    1    2    3    5    8   13   21   34   55   89  144  233  377  610  987 1597

Loops können mit break unterbrochen werden - das gilt nicht nur für repeat, sondern auch für die anderen beiden Formen von Loops. Hier wurde eine if-Bedingung in den Loop geschachtelt. In jedem einzelnen Durchlauf des Loops wird geprüft, ob die Bedingung erfüllt ist, und die Durchführung wird beendet (break), sobald dies der Fall ist.

Ergänzen Sie print(fibo) vor der if-Abfrage, und schauen Sie sich das Ergebnis an. Dies zeigt Ihnen gewissermaßen das “Innenleben” Ihres Loops. Sie sehen so genauer, was in jedem Schritt des Loops passiert, und können oftmals leichter nachvollziehen, wodurch beispielsweise Fehler entstehen.

An den folgenden Laufzeiten sehen wir, dass Loops tatsächlich deutlich langsamer sind, als die vektorwertige Alternative.

# simuliere 1000 Beobachtungen und bestimme den Mittelwert 
X <- rnorm(10^3, mean = 1, sd = 2)
m <- mean(X)

t1 <- Sys.time() # speichere die Startzeit
X_c <- X
for(x in X)
{
  X_c[i] <- x - m
}
Sys.time() - t1 # bestimmte die Laufzeit durch aktuelle Zeit minus Startzeit

## Time difference of 0.01244545 secs

t2 <- Sys.time() # speichere die Startzeit
X_c <- X - m
Sys.time() - t2 # bestimmte die Laufzeit durch aktuelle Zeit minus Startzeit

## Time difference of 0.001213074 secs

Loops sind in diesem Beispiel fast um den Faktor 10 langsamer (zumindest, wenn Sie den Code für 10^6 Beobachtungen durchführen, die Maske, die für diese Website genutzt wird, ist deutlich langsamer, weswegen hier auf 10^3 ausgewichen wurde…). Es gibt jedoch viele Anwendungsgebiete, wo Loops das Mittel der Wahl sind!

apply-Funktionen

Auch apply und seine Varianten können genutzt werden, um bspw. einen for-Loop auszudrücken. Diese Funktion verkürzt die Schreibweise und kann manchmal auch die Laufzeit verkürzen, insbesondere wenn bspw. das pbapply-Paket verwendet wird, welches einfaches Parallelisieren erlaubt.

A <- data.frame("a" = c(2,3,4), "b" = c(1,1,1))
apply(A, 2, mean) # Mittelwert über Spalten/Variablen

## a b 
## 3 1

colMeans(A)

## a b 
## 3 1

apply(A, 1, mean) # Mittelwert über Zeilen/Personen/Beobachtungen

## [1] 1.5 2.0 2.5

rowMeans(A)

## [1] 1.5 2.0 2.5

apply(A, 2, sd) # Standardabweichung über Spalten/Variable

## a b 
## 1 0

In der Sitzung zu Simulationsstudien und Poweranalysen aus dem vergangenen Semester hatten wir empirisch die Power und den $\alpha$-Fehler des $t$-Tests sowie des Korrelationstest untersucht. Dabei hatten wir replicate verwendet. Bspw. hatten wir mit folgendem Code den $p$-Wert des $t$-Tests unter der $H_0$ Hypothese untersucht:

N <- 20
set.seed(1234)
replicate(n = 10, expr = {X <- rnorm(N)
                          Y <- rnorm(N)
                          ttestH0 <- t.test(X, Y, var.equal = TRUE)
                          ttestH0$p.value})

##  [1] 0.26352442 0.03081077 0.21285027 0.27429670 0.53201656 0.79232864 0.93976306 0.43862992 0.96766599
## [10] 0.68865560

Wenn wir nun genauer hinschauen, dann sehen wir, dass der Block

{X <- rnorm(N)
 Y <- rnorm(N)
 ttestH0 <- t.test(X, Y, var.equal = TRUE)
 ttestH0$p.value}

im Grunde nichts weiter darstellt, als das Innere einer Funktion. replicate ist im Grunde nichts anderes als eine bestimmte for-Schleife, nämlich eine for-Schleife, in welcher das Argument nicht genutzt wird! Wir schreiben das Ganze mal mittels einer Funktion:

mySim <- function(N)
{
  X <- rnorm(N)
  Y <- rnorm(N)
  ttestH0 <- t.test(X, Y, var.equal = TRUE)
  return(ttestH0$p.value)
}
set.seed(1234)
replicate(n = 10, expr = mySim(N = 20))

##  [1] 0.26352442 0.03081077 0.21285027 0.27429670 0.53201656 0.79232864 0.93976306 0.43862992 0.96766599
## [10] 0.68865560

In der Sitzung zu Simulationsstudien und Poweranalysen hatten wir außerdem den empirischen $t$-Wert untersucht. Diesen können wir nun ganz leicht mit aufnehmen.

mySim2 <- function(N)
{
  X <- rnorm(N)
  Y <- rnorm(N)
  ttestH0 <- t.test(X, Y, var.equal = TRUE)
  return(c("p" = ttestH0$p.value, "t" = ttestH0$statistic))
}
set.seed(1234)
replicate(n = 10, expr = mySim2(N = 20))

##          [,1]        [,2]      [,3]       [,4]      [,5]      [,6]      [,7]       [,8]        [,9]
## p   0.2635244  0.03081077 0.2128503  0.2742967 0.5320166 0.7923286 0.9397631  0.4386299  0.96766599
## t.t 1.1349024 -2.24295556 1.2670437 -1.1092419 0.6306927 0.2651479 0.0760693 -0.7827414 -0.04080374
##         [,10]
## p   0.6886556
## t.t 0.4037557

Wir sehen, dass die p-Werte und die t-Werte nun gleichzeitig ausgegeben werden und zwar in zwei Zeilen untereinander, da wir den Output als Vektor gewählt haben! In diesem Semester hatten wir uns bisher mit der Regressionsanalyse beschäftigt. Aus diesem Grund wollen wir an dieser Stelle noch kurz anschneiden, wie eine Simulationsstudie für eine Regression durchgeführt werden könnte. Zunächst brauchen wir dazu Prädiktoren. Mit Hilfe der rmvnorm Funktion aus dem mvtnorm-Paket lassen sich leicht multivariat-normalverteilte Zufallsvariablen simulieren, deren Mittelwerte und Kovarianz bekannt ist:

S <- matrix(c(1, .7, .7, 2), 2, 2) # Populationskovarianzmatrix
S

##      [,1] [,2]
## [1,]  1.0  0.7
## [2,]  0.7  2.0

# install.packages("mvtnorm")
library(mvtnorm)
set.seed(1234)
X <- rmvnorm(n = 10^3, mean = c(2, 3), sigma = S)
colMeans(X)

## [1] 1.997926 2.980730

cov(X)

##          [,1]     [,2]
## [1,] 1.046158 0.719314
## [2,] 0.719314 1.871493

Für größeres n landen wir näher bei den Populationswerten:

$$ \mathbb{E}[X] = \begin{pmatrix} 2 \ 3 \end{pmatrix},\quad \mathbb{C}ov[X]=\begin{pmatrix} 1 & .7 \ .7 & 2 \end{pmatrix}.$$

Für eine Regressionsanalyse brauchen wir jetzt nur noch ein Residuum, sowie die $\beta$-Gewichte, um die Variable $Y$ zu definieren. Angenommen wir wollen folgendes Populationsmodell untersuchen:

$$Y_i = 0.3 + 0.5\cdot X_{1i} + 0.3\cdot X_{2i} + \varepsilon_i$$

wobei $\varepsilon_i$ eine Residualstandardabweichung von 1.3 haben soll:

eps <- rnorm(10^3, sd = 1.3)
X1 <- X[,1]
X2 <- X[,2]
Y <- 0.3 + 0.5*X1 + 0.3*X2 + eps
df <- data.frame("X1" = X1, "X2" = X2, "Y" = Y)

Dann können wir nun leicht eine Regressionsanalyse durchführen:

reg <- lm(Y ~ 1 + X1 + X2, data = df)
coef(reg) # Koeffizienten abgreifen

## (Intercept)          X1          X2 
##   0.4480455   0.5145347   0.2532168

Wir sehen, dass die Koeffizienten recht nah an den “wahren” Werten liegen. Verpacken wir das Ganze in eine Funktion, so können wir den Bias der Schätzung untersuchen. Der Bias ist die durchschnittliche Abweichung der Schätzung vom wahren Wert. Ein Bias von 0 ist somit erstrebenswert!

myRegSim <- function(N)
{
  S <- matrix(c(1, .7, .7, 2), 2, 2) # Populationskovarianzmatrix
  X <- rmvnorm(n = N, mean = c(2, 3), sigma = S)
  eps <- rnorm(N, sd = 1.3)
  X1 <- X[,1]
  X2 <- X[,2]
  Y <- 0.3 + 0.5*X1 + 0.3*X2 + eps
  df <- data.frame("X1" = X1, "X2" = X2, "Y" = Y)
  reg <- lm(Y ~ 1 + X1 + X2, data = df)
  coef(reg) # Koeffizienten abgreifen
  return(coef(reg))
}
set.seed(1234)
replicate(n = 10, expr = myRegSim(N = 10^3))

##                  [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]      [,9]
## (Intercept) 0.4480455 0.4645526 0.0959823 0.4036081 0.3621404 0.4454766 0.2000509 0.2704179 0.2343473
## X1          0.5145347 0.4229675 0.5794930 0.6167309 0.4935631 0.4068533 0.5341167 0.4260486 0.5119222
## X2          0.2532168 0.3155510 0.3264700 0.2033327 0.2789048 0.3184694 0.3038769 0.3344544 0.3198885
##                 [,10]
## (Intercept) 0.4723725
## X1          0.4416930
## X2          0.3067544

Speichern wir das Ganze ab, transponieren es und bilden colMeans, so erhalten wir eine Schätzung für die durchschnittliche Schätzung unseres Experiments (das wir insgesamt 10 Mal unter identischen Voraussetzungen durchführen konnten):

set.seed(1234)
mySimErg <- t(replicate(n = 10, expr = myRegSim(N = 10^3)))
colMeans(mySimErg)

## (Intercept)          X1          X2 
##   0.3396994   0.4947923   0.2960919

Selbst bei nur 10 Wiederholungen und einer Stichprobengröße von 1000 ist der Bias schon sehr gering (zum Vergleich, die wahren Werte waren 0.3, 0.5, 0.3). Der Bias wird nun so bestimmt:

$$\hat{\theta} - \theta_0,$$

wobei $\hat{\theta}$ der durchschnittliche Koeffizient und $\theta_0$ der wahre Koeffizient ist. In den Übungen werden wir uns mit diesem Sachverhalt noch etwas genauer beschäftigen und die neuen Erkenntnisse der Regressionsanalyse mit einfachen Simulationen untermauern. Weitere interessante Größen einer Simulationsstudie sind die Power und der $\alpha$-Fehler, sowie der mittlere Standardfehler und der Vergleich zwischen mittlerem Standardfehler ($\bar{SE}(\hat{\theta})$ mittlerer geschätzter Streuung der Schätzungen anhand der Daten) und der tatsächlichen Streuung der Schätzungen ($SD(\hat{\theta})$). Diese Werte können wir untersuchen, indem wir die summary auf das Regressionsobjekt anwenden und anschließend wieder mit coef die Parameterschätzer Estimate, die Standardfehler Std.Error, sowie den $t$-Wert und den zugehörigen $p$-Wert bestimmen. Diese Werte sind entscheidend für Signifikanzentscheidungen und somit für Power und $\alpha$-Fehler.

Zum Abschluss noch ein Gedankenexperiment: Wenn wir immer wieder Daten simulieren, dann erhalten wir von Mal zu Mal unterschiedliche Parameterschätzer. Die Estimates streuen also. Diese Streuung beschreibt die Unsicherheit, die beim Schätzen in einem “endlichen” Sample entsteht. Sie nennt sich bezogen auf eine Monte-Carlo-Simulationsstudie (die wir gerade durchgeführt hatten) MCSD (für Monte-Carlo-Standardabweichung). Der Standardfehler, den wir mit jedem Mal Schätzen bekommen, soll nun eine Schätzung für diese wahre Streuung der Parameterschätzungen sein. Aus diesem Grund wird in Simulationsstudien häufig der durchschnittliche Standardfehler (MCSE) mit der MCSD verglichen. Eine gute Schätzung des SEs für MCSD ist also entscheidend dafür, dass der statistische Test vertrauenserweckende Ergebnisse liefert!