In diesem Block werden wir verschiedene Arten von Loops (Schleifen) kennenlernen und lernen, vertieft mit Funktionen zu arbeiten. Dieses Wissen wollen wir dann nutzen, um nochmals Power- bzw. Simulationsanalysen durchzuführen, welche wir in der Sitzung zu Simulationsstudien und Poweranalysen bereits kennengelernt haben. Vorab beschäftigen wir uns noch mit einigen Grundlagen zum Thema logische Abfragen. Show
Logische Abfragen und Bedingungen: if und elseIm Prozess der Datenaufbereitung und -auswertung kommt man häufig an den Punkt, an dem ein bestimmter Befehl nur unter bestimmten Bedingungen ausgeführt werden soll, oder in dem abhängig von einer Bedingung unterschiedliche Aktionen ausgeführt werden sollen. Dabei bezieht sich die Bedingung auf einen Wert in einer bestimmten Variable, der sich zwischen den Versuchspersonen unterscheidet. Dafür können wir sogenannte Wenn-Dann-Bedingungen, oder auch if-Abfragen nutzen, in denen wir definieren, unter welchen Bedingungen ein folgender Befehl ausgeführt werden soll. >Beispiel: In einer neuen Variable Med wollen wir für alle Versuchspersonen eine 1 vergeben, die in der Variable “Dosis” einen gültigen Wert haben, und eine 0 vergeben für alle Personen, die in der Variable “Dosis” ein NA haben. Im letzten Semester hatten wir das über die
Bedingungen in eckigen Klammern angewandt auf einen if-AbfragenEinfache if-AbfrageWie in eigentlich allen Programmiersprachen werden Wenn-Dann-Bedingungen auch in Die Konsequenz wird nur ausgeführt, wenn die Bedingung das Ergebnis
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
In diesem Fall stimmt die logische Abfrage ( Wenn jedoch
Wenn Sie diesen Code wieder Schritt für Schritt ausführen (z.b. lediglich die Abfrage in den runden Klammern), können Sie sich so die Zusammenhänge innerhalb der if-Abfrage verdeutlichen. if-Abfrage mit mehreren MöglichkeitenOft 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
Hier erhalten
wir die Antwort, ja, Monica ist eine Figur aus der Serie. Der Ausdruck Wenn wir die gleiche Abfrage auf eine andere Person anwenden, trifft die Bedingung nicht zu, und der Befehl wird nicht ausgeführt.
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
An diesem Beispiel wird der Unterschied der beiden Befehle deutlich. Um das Ganze nochmals zu untermauern:
Mit diesen beiden Befehlen können wir leicht prüfen, ob alle Elemente oder mindestens einer eine Eigenschaft erfüllt. Abgleich mit einem DatumEs sind beispielsweise auch logische Abfragen mit Zeitpunkten und Daten möglich. Zum Beispiel
können wir mit dem Befehl
Verknüpfung logischer AbfragenWie im letzten Semester bereits
besprochen, können logische Bedingungen mit
Durch die logische Verknüpfung mit 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-AbfragenHä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. Um zwischen zwei alternativen Befehlen auszuwählen, ergänzen wir das
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
Das
gilt sowohl für else if-BedingungenHä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 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.
Wir versuchen nachzuvollziehen, was in dieser verschachtelten
Wenn dem so ist, wird der nächste Block ausgeführt:
Dieser fragt ab, ob heute Montag ist ( Ist heute kein Wochentag, dann wird direkt die
Dieses Spiel der Verschachtelung lässt sich beliebig erweitern. Natürlich sind Wenn-Dann-Abfragen eigentlich hauptsächlich dann nützlich, wenn der Code für verschiedene Daten, Objekte oder Funktionen mehrfach genutzt werden soll und man nicht in jedem Einzelfall schon vorher weiß, welche Inhalte die Objekte haben, mit denen man arbeitet. Ein einfaches Beispiel mit einer zufällig gezogenen Zahl könnte so aussehen:
Wenn Sie diesen Code mehrfach ausführen, bekommen Sie immer wieder unterschiedliche Paare an Funktion ifelseWenn 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
So wird die gleiche
Loops (Schleifen)Beim Programmieren kommt es häufig vor, dass der gleiche Befehl mehrfach angewandt werden muss. Loops (oder Schleifen) bieten die Möglichkeit, den gleichen
Gerade in Kombination mit
|
Variable | Adjektiv | Richtung | Dimension |
---|---|---|---|
stim1
| zufrieden | positiv | Gut vs. Schlecht |
stim2
| ausgeruht | positiv | Wach vs. Müde |
stim3
| ruhelos | negativ | Ruhig vs. Unruhig |
stim4
| schlecht | negativ | Gut vs. Schlecht |
stim5
| schlapp | negativ | Wach vs. Müde |
stim6
| gelassen | positiv | Ruhig vs. Unruhig |
stim7
| müde | negativ | Wach vs. Müde |
stim8
| gut | positiv | Gut vs. Schlecht |
stim9
| unruhig | negativ | Ruhig vs. Unruhig |
stim10
| munter | positiv | Wach vs. Müde |
stim11
| unwohl | negativ | Gut vs. Schlecht |
stim12
| entspannt | positiv | Ruhig vs. Unruhig |
In der Spalte Dimension sehen wir, dass die Items 3 verschiedene Dimensionen abbilden: Gut vs. Schlecht, Wach vs. Müde und Ruhig vs. Unruhig. Die Items sind dabei unterschiedlich gepolt - die Adjektive “ausgeruht” und “schlapp” erfasst beide die Dimension Wach vs. Müde, jedoch in unterschiedlicher Ausrichtung. Um die drei Skalenwerte berechnen zu können müssen die jeweils “negativen” Adjektive ins Positive umgepolt werden. Hierzu gibt es zum Beispiel folgende zwei Möglichkeiten, die wir bereits aus dem vergangenen Semester kennen. Zum Einen können wir bei den entsprechenden Items die Skalenwerte ersetzen:
mdbf$stim4_r[mdbf$stim4 == 1] <- 4
mdbf$stim4_r[mdbf$stim4 == 2] <- 3
mdbf$stim4_r[mdbf$stim4 == 3] <- 2
mdbf$stim4_r[mdbf$stim4 == 4] <- 1
Oder wir können das Vorgehen verkürzen, indem wir die folgende Berechnungsweise anwenden:
mdbf$stim4_r <- -1 * (mdbf$stim4 - 5)
Aber trotz der Verkürzung haben wir nun erst ein einziges Item umcodiert. Mit Hilfe von Loops können wir uns die Arbeit ersparen, diesen Abschnitt für jedes negative
Adjektiv schreiben zu müssen. Wir erinnern uns: Für den for
-Loop müssen wir wissen, für welche Fälle ein Skript durchgeführt werden muss. Für die Umcodierung der Items speichern wir also alle negativen Items in einem Vektor neg
(Häufig werden auch die Spaltennummern verwendet):
# Kopie des Datensatzes erstellen, um Datenverlust vorzubeugen
mdbf_r <- mdbf
# Vektor der negativen Items
neg <- c("stim3", "stim4", "stim5", "stim7", "stim9", "stim11")
In neg
wird kodiert, welche Items negativ formuliert sind, und in die Umcodierung einbezogen werden sollen. Danach wenden wir die oben bereits gezeigte Formel erneut an, hier jedoch
nacheinander auf jedes der Elemente, die in neg
gespeichert sind. Dabei nimmt i nacheinander die Namen der Variablen an stim3, stim4, stim5, stim7,… an. Nun können wir mit mdbf_r[, i]
(bzw. mdbf[, i]
) die richtige Spalte im Datensatz anwählen (mbdf_r[[i]]
geht entsprechend auch, da es sich bei mdbf_r
um einen data.frame
handelt), also im ersten Schritt mit mdbf_r[, "stim3"]
(bzw. mdbf[, "stim3"]
) die dritte Variable stim3, also das dritte Item. Der Platzhalter i
iteriert also durch die Elemente von
neg
. Damit bei mehrmaligem durchführen des Skriptes nicht immer hin- und herkodiert wird, speichern wir die Berechnung mit Hilfe von mdbf[, i]
in mdbf_r[, i]
ab.
for (i in neg) {
mdbf_r[, i] <- -1 * (mdbf[, i] - 5)
}
Zur Prüfung des Erfolges berechnen wir die Korrelation des Items stim3
im originalen Datensatz und im umcodierten Zustand.
cor(mdbf[, "stim3"], mdbf_r[, "stim3"])
## [1] -1
Um Ihr Verständnis zu überprüfen, versuchen Sie, in einer neuen Kopie des Datensatzes jetzt stattdessen alle positiven Items umzucodieren!
Zudem können for
-Loops
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.
Weitere Loops
Weitere häufig verwendete Loops
sind die while
und die repeat
-Loop: In while
-Loops wird der Code so lange ausgeführt, bis eine vorab definierte Bedingung erfüllt ist. 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. Eine Loop (auch for
oder while
) können mit dem break
-Befehl beendet werden. Für mehre Informationen zu while
, repeat
und break
siehe Appendix A.
Anmerkung: Generell sollten Loops in R
nur genutzt werden, wenn keine vektorbasierte Alternative zur Verfügung steht. Zum Beispiel:
um eine Variable zu zentrieren sollte nicht ein Loop genutzt werden, der von jedem Element des Vektors den Mittelwert abzieht (bei Interesse an einem Laufzeitvergleich siehe Appendix B. Stattdessen ist R
in der Lage den Mittelwert direkt von jedem Element des Vektors abzuziehen (elementeweise Anwendung) - diese Umsetzung ist also direkt Vektor-basiert und in R
(beinahe ausnahmslos) die
schnellere und effizientere Variante.
Funktionen
Sie haben bereits gelernt, dass (fast) alle Aktionen, die in R
ausgeführt werden, sich sogenannte Funktionen zunutze machen. Hier wollen wir noch einen Schritt weiter gehen, und lernen, wie Sie selbst Funktionen schreiben können. Funktionen, die in R
angewendet werden können, sind ebenfalls Objekte. Dadurch können eigene Funktionen wie andere Objekte auch angelegt werden - dazu
müssen sie lediglich mit der function
-Funktion erstellt werden. Im Allgemeinen sieht das wie folgt aus:
eigene_funktion <- function(argument1, argument2, ...) {
# Durchgeführte Operationen
}
Der Name der erstellen Funktion steht hier ganz am Anfang. function
ist die Funktion, die dafür zuständig ist, neue Funktionen zu definieren. In den runden Klammern dahinter müssen Sie angeben, welche Argumente Ihre Funktion annehmen soll. Auf diese Argumente können Sie in der Beschreibung der Operationen zugreifen. In geschweiften Klammern geben Sie als nächstes an,
welche Operationen mit den genannten Argumenten durchgeführt werden sollen. Als Argumente können beliebig viele Einstellungen für die Funktion definiert werden, auf die dann in der Funktion Bezug genommen wird. Wichtig ist dabei, dass Funktionen keinen generellen Zugriff auf den Workspace haben, sondern alle Objekte, die sie benötigen, durch die Argumente an sie weitergegeben werden müssen. Auch innerhalb einer Funktion definierte Objekte können nicht außerhalb der Funktion genutzt werden!
Beispiel Varianzfunktion
In R
wird mit der var
-Funktion die Schätzung für die Populationsvarianz \(\widehat{\sigma}^2\) und nicht die empirische Varianz \(s\) bestimmt. Wir könnten also eine eigene Funktion anlegen, die die empirische Varianz schätzt. Dafür können wir die Formel zur Varianzberechnung einfach in R
-Code übersetzen:
\[s^2 = \frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n}\]
Als
R
-Code würde wir also zunächst die einzelnen Elemente definieren, und dann nach dem Vorbild der Formel die Varianz berechnen.
x <- mdbf[, 1]
n <- length(x)
s2 <- sum((x - mean(x))^2) / n
s2
## [1] 0.482299
Dieser Code funktioniert allerdings nur für eine einzige Variable mit dem Namen x
und wir müssten den Code für jede einzelne Anwendung wiederholen. Um das abzukürzen, können wir eine eigene, wiederverwendbare Funktion anlegen. Dafür nutzen wir wie oben beschrieben die Funktion function
. Unsere neue Funktion soll empVar
heißen,
und erhält nur ein einziges Argument x
. In den geschweiften Klammern definieren wir, wie die Berechnung funktionieren soll.
empVar <- function(x) {
n <- length(x)
s2 <- sum((x - mean(x))^2)/n
}
empVar(mdbf[, 1])
Nun erhalten wir jedoch kein Ergebnis, wenn wir diese Funktion auf empVar(mdbf[, 1])
anwenden. Dafür müssen wir zusätzlich mit return
definieren, was der Ausgabewert der Funktion sein soll. In diesem Fall wird das Ergebnis der Berechnung ausgegeben. Wenn kein return
-Wert definiert wird, gibt die Funktion bei der Anwendung kein Ergebnis in die Konsole
aus. Wir haben auch keinen Zugriff auf das Objekt s2. Eine Funktion ohne return
wird zwar ausgeführt, man hat aber keinen Zugriff auf das Ergebnis, weil alle innerhalb der Funktion angelegten Objekte entfernt werden, sobald die Durchführung der Funktion abgeschlossen ist. Funktionen sollten also prinzipiell mit return
Ergebnisse nach außen geben:
empVar <- function(x) {
n <- length(x)
s2 <- sum((x - mean(x))^2)/n
return(s2)
}
Diese Funktion kann jetzt auf jede beliebige Variable angewendet werden:
empVar(mdbf[, 1])
## [1] 0.482299
empVar(mdbf[, 2])
## [1] 0.7213661
Das Einzige, was diese
Funktion von in R
implementierten Paketen unterscheidet ist, dass sie explizit im Workspace bzw. Environment angezeigt wird. Dies können Sie mit dem ls()
-Befehl ausprobieren. Außerdem sehen Sie die Funktion in R
-Studio oben rechts in dem Fenster, in dem auch die Daten etc. dargestellt werden unter Functions.
Weil beim Durchführen von Funktionen als erstes der Workspace nach definierten Funktionen durchsucht wird, sollten Funktionen möglichst einzigartig
benannt werden, weil sonst nicht mehr (so leicht) auf die R
-internen Funktionen zugegriffen werden kann.
Zusätzlich sollte beachtet werden, dass return
nur ein einziges Argument entgegennimmt: Funktionen in R
können also nur ein einziges Objekt als Ergebnis liefern. Wenn mehrere Ergebnisse ausgegeben werden sollen, müssen diese vorher innerhalb der Funktion zu einem Objekt (meistens einer Liste) zusammengefasst werden. Eine Liste kann benannt werden, indem zunächst die
Namen geschrieben werden und diesen dann via =
etwas zugeordnet wird:
empVar <- function(x) {
n <- length(x)
s2 <- sum((x - mean(x))^2)/n
out <- list("s2" = s2, "n" = n)
return(out)
}
empVar(mdbf[, 2])
## $s2
## [1] 0.7213661
##
## $n
## [1] 98
Funktionen können eine beliebige Anzahl von Argumenten entgegennehmen, aber nur ein einziges Objekt als Ergebnis liefern (wobei dieses Objekt natürlich mehrere Elemente haben kann, bspw. hat eine Vektor ja mehrere Einträge oder eine Liste kann auch aus mehreren Vektoren, Matrizen oder anderem bestehen). Um eine gemeinsame Funktion für beide Formen der Varianz zu haben, könnten wir die Anzahl der Argumente erweitern.
Vari <- function(x, empirical) {
n <- length(x)
if (empirical) {
s2 <- sum((x - mean(x))^2)/n
} else {
s2 <- sum((x - mean(x))^2)/(n-1)
}
return(s2)
}
Das Argument namens empirical
kann in dieser Funktion genutzt werden, um zu entscheiden, welche Varianz-Formel angewandt werden soll. In diesem Fall wird also eine Einstellung für empirical
benötigt, die dann von if
als TRUE
oder FALSE
bewertet werden kann:
Vari(mdbf[, 2], TRUE)
## [1] 0.7213661
Vari(mdbf[, 2], FALSE)
## [1] 0.7288029
Wenn wir die Einstellung vergessen, wird - wie bei allen anderen R
Funktionen auch - ein Fehler produziert. Probieren Sie dies aus, und beachten Sie den
Wortlaut der Fehlermeldung. Jetzt sollten Sie verstehen, was diese aussagt und wie der Fehler behoben werden kann!
Vari(mdbf[, 2])
## Error in Vari(mdbf[, 2]) :
## argument "empirical" is missing, with no default
Die Fehlermeldung beinhaltet die Worte with no default. Dies impliziert, dass eine Voreinstellung für das Argument gesetzt werden könnte. Dann müssen Nutzer:innen das Argument nicht mehr zwingend angeben. Wenn Voreinstellungen für Argumente festgelegt werden sollen, erreichen wir das, indem in der runden Klammer direkt der default-Wert für ein Argument mit angegeben wird.
Vari <- function(x, empirical = TRUE) {
n <- length(x)
if (empirical) {
s2 <- sum((x - mean(x))^2)/n
} else {
s2 <- sum((x - mean(x))^2)/(n-1)
}
return(s2)
}
Vari(mdbf[, 2])
## [1] 0.7213661
Solange jetzt nicht explizit etwas bei der Anwendung der Funktion für empirical
deklariert wird, wird von der Voreinstellung TRUE
ausgegangen.
Außerdem gibt es noch Funktionen wie replicate
, apply
, lapply
, sapply
, etc. die die Grenzen zwischen Schleifen und Funktionen etwas verwaschen. Siehe auch Appendix A.
Anwendung: Simulationsstudien und Poweranalysen
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
## [7] 0.93976306 0.43862992 0.96766599 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
## [7] 0.93976306 0.43862992 0.96766599 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]
## p 0.2635244 0.03081077 0.2128503 0.2742967 0.5320166 0.7923286 0.9397631
## t.t 1.1349024 -2.24295556 1.2670437 -1.1092419 0.6306927 0.2651479 0.0760693
## [,8] [,9] [,10]
## p 0.4386299 0.96766599 0.6886556
## t.t -0.7827414 -0.04080374 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]
## (Intercept) 0.4480455 0.4645526 0.0959823 0.4036081 0.3621404 0.4454766
## X1 0.5145347 0.4229675 0.5794930 0.6167309 0.4935631 0.4068533
## X2 0.2532168 0.3155510 0.3264700 0.2033327 0.2789048 0.3184694
## [,7] [,8] [,9] [,10]
## (Intercept) 0.2000509 0.2704179 0.2343473 0.4723725
## X1 0.5341167 0.4260486 0.5119222 0.4416930
## X2 0.3038769 0.3344544 0.3198885 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!
Appendix A
while
und repeat
-Loops & weitere Loop-Funktionenwhile-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" "Kopf" "Kopf" "Zahl" "Zahl" "Zahl" "Kopf" "Kopf" "Zahl" "Kopf"
## [11] "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
## [16] 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.
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