R Befehle (functions)

Bekanntlich ist alles in R ein Objekt, auch R Befehle (functions), sind davon keine Ausnahme. Eine sehr gute Einführung bieten Reto Stauffer, Thorsten Simon, Achim Zeileis, siehe https://eeecon.uibk.ac.at/~discdown/rprogramming/functions.html

Grundlagen

Alles, was in R passiert, geschieht durch einen Funktionen-Aufruf (call).

Funktionen haben in den meisten Fällen Inputs und Outputs, wenn wir z.B. der Funktion mean() die Zahlen 1 bis 5 übergeben erhalten wir den Mittelwert 3 zurück

mean(1:5)

## [1] 3

In den runden Klammern, die unmittelbar nach dem Funktionsnamen (d.h. ohne Leerzeichen!) folgen müssen, kann eine Argumentenliste übergeben werden, im obigen Fall ein Vektor mi den Zahlen 1 bis 5.

Nicht alle Funktionen müssen einen Input haben, z.B. können wir mit q() (für quit) R beenden. Wenn wir einen Funktionsnamen ohne Klammern eingeben erhalten wir den Code (bzw. Inhalt) der Funktion (oder zumindest einen Teil davon)

## function (save = "default", status = 0, runLast = TRUE) 
## .Internal(quit(save, status, runLast))
## <bytecode: 0x0000000013928f58>
## <environment: namespace:base>

Dies zeigt uns, dass die Funktion q() drei Argumente hat, save, status und runLast. Das erste Argument save erlaubt zu bestimmen, ob der workspace gespeichert werden soll. Wie wir sehen ist diesem Argument als Weglass- (default) Wert (ein character string) "default" zugeordnet (d.h. die Konfigurations-Einstellungen werden übernommen); alternativ hätten wir "no", "yes" oder "ask" wählen können.

Welche Argumente und Optionen zur Verfügung stehen erfährt man z.B. mit der help(q) Funktion, oder einfacher, ?q.

Mit

q(save = "no")  ## oder einfacher 
q("no")

können wir bestimmen, dass der workspace nicht gespeichert wird. Die zweite Variante q("no") funktioniert, weil save das erste Argument der Liste ist, und diesem wird der String "no"zugewiesen. Die einzelnen Argumente werden durch einen Beistrich getrennt.

Das zweite Argument status gibt eine numerischen Fehlercode zurück (0 bedeutet erfolgreiche Ausführung). Das dritte Argument runlast kann entweder TRUE oder FALSE sein und gibt an, ob vor dem Beenden noch bestimmte Operationen durchgeführt werden sollten.

Mit

q(runLast = FALSE) ## oder 
q( , , FALSE)

können wir angeben, dass diese Operationen nicht ausgeführt werden sollen. Beachte, dass bei der zweiten Variante die Positon an dritter Stelle (durch die Beistriche definiert) angegeben werden muss. Dies ist nicht nur unübersichtlich, sondern auch sehr fehleranfällig, deshalb sollte immer die erste Variante mit Angabe der Namen der Argumente gewählt werden!

Einer der vielen Vorteile von RStudio ist, dass man durch Drücken der Tab-Taste in den geschlossenen Klammern der Funktion eine Liste der verfügbaren Argumente mit nützlichen Hinweisen vorgeschlagen bekommt.

Achtung: Innerhalb von Argumentenlisten von Funktionen erfolgt die Zuweisung mit mit dem = Zeichen, nicht (!) mit einem Zuweisungspfeil <-. Warum? Würden wir z.B. q(runlast <- FALSE) eingeben, würde R ein neues Objekt mit dem Namen runlast anlegen, und diesem neuen Objekt FALSE zuordnen. Da dies an erster Stelle steht wäre dies gleichbedeutend mit q(save = FALSE), vermutlich nicht das, was wir wollen!

Eigene Funktionen

Eine große Stärke von R ist, dass man mit function() auch sehr einfach eigene Funktionen (Befehle) schreiben kann, z.B.

hw <- function() {return(print("Hallo Welt!"))} 
hw()

## [1] "Hallo Welt!"

Man beachte, dass für die Bildschirmausgabe innerhalb von Funktionen der print() Befehl erforderlich ist.

Schon etwas komplizierter, eine (ziemlich sinnlose) Funktion, die zu jeder eingegeben Zahl 3 addiert

add3 <- function(x) { 
  z <- x + 3 
  return(z) 
}

und wenn wir diese neue Funktion ausführen

add3(2)
## [1] 5

Wir könnten auch noch mit einer if Abfrage überprüfen, ob tatsächlich ein numerischer Wert übergeben wurde, und anderenfalls eine Fehlermeldung (not a number) zurückgegeben.

add3 <- function(x) { 
  if (is.numeric(x)) {
    z <- x + 3 
  } else { 
    z <- "not a number" 
  } 
  return(z) 
}

und

add3(5) 
## [1] 8
add3("Apfel") 
## [1] "not a number"

Man kann auch in eine existierende R-Funktion eine eigene Funktion einfügen. Zum Beispiel die Funktion sapply(list, FUN), womit auf jedes Element einer Liste (meist eines data.frame) eine Funktion FUN angewandt wird. Diese Funktion ersetzt im wesentlichen eine (meist viel langsamere) Schleife.

Diese geschachtelten Funktionen können wir z.B. verwenden, um zu für alle Variablen eines data.frame die Anzahl der gültigen Beobachtungen zu zählen.

Wir legen zuerst einen data.frame mit drei Variablen an, wobei die zweite und dritte Variable (y und z) jeweils missings (NA) enthalten

df <- data.frame( x = c(1, 2, 3, 4, 5),
                  y = c(1, NA, 3, NA, 5), 
                  z = c(1, NA, NA, NA, NA) ) 
summary(df)

##        x           y           z    
##  Min.   :1   Min.   :1   Min.   :1  
##  1st Qu.:2   1st Qu.:2   1st Qu.:1  
##  Median :3   Median :3   Median :1  
##  Mean   :3   Mean   :3   Mean   :1  
##  3rd Qu.:4   3rd Qu.:4   3rd Qu.:1  
##  Max.   :5   Max.   :5   Max.   :1  
##              NA's   :2   NA's   :4

Der summary() Befehl sagt uns zwar, wie viele Werte fehlen (also NA für not available sind), aber nicht unmittelbar, wieviele gültige Einträge vorliegen.

Mit der Funktion

is.na(df$y)

## [1] FALSE  TRUE FALSE  TRUE FALSE

erhalten wir einen logischen Vektor zurück, der TRUE für missings ist; mit der Negation !is.na(df$y) sind gültige Einträge TRUE, und diese können wir mit sum(!is.na(df$y)) zählen.

sum(!is.na(df$y))

## [1] 3

Wenn dies für alle Variablen eines data.frame geschehen soll können wir die oben erwähnte `sapply()’ Funktion verwenden

sapply(df, function(x) {sum(!is.na(x))})

## x y z 
## 5 3 1

Der `lm()` Befehl

Da wir die lm() Funktion (für linear model) im Folgenden sehr häufig benötigen werden, einige Anmerkungen dazu. Er hat die Form eq <- lm(formula, data, ...), wobei formula ein Ausdruck der Art y ~ x1 + x2 ist, und das Symbol ~ wird gelesen als ‘wird erklärt durch’. Die drei Punkte ... stehen für weitere Argumente.

In formula Objekten hat das + nicht die übliche arithmetische Bedeutung, sondern stellt einen linearen Zusammenhang dar. Ähnlich haben auch die anderen Rechenzeichen innerhalb von formula Objekten eine andere Bedeutung, die v.a. für die Eingabe von Interaktionseffekten nützlich ist.

Zum Beispiel erzeugt y ~ x1 * x2 das gleiche wie y ~ x1 + x2 + I(x1 * x2), wobei im zweiten Ausdruck die I() Funktion (für inhibit) benötigt wird, damit wir das arithmetische Multiplikationszeichen in der formula Umgebung ‘geschützt’ und richtig interpretiert wird, z.B.

set.seed(12345)
x1 <- rnorm(n = 10)
x2 <- rnorm(n = 10)
y  <- 5*x1 - 2*x2 + rnorm(n = 10)
lm(y ~ x1 * x2)

## 
## Call:
## lm(formula = y ~ x1 * x2)
## 
## Coefficients:
## (Intercept)           x1           x2        x1:x2  
##     -0.1710       4.8686      -1.2486      -0.1074

lm(y ~ x1 + x2 + I(x1 * x2))

## 
## Call:
## lm(formula = y ~ x1 + x2 + I(x1 * x2))
## 
## Coefficients:
## (Intercept)           x1           x2   I(x1 * x2)  
##     -0.1710       4.8686      -1.2486      -0.1074

Ebenso haben in formula Objekten die Zeichen /, ^, :, * eine spezielle Bedeutung, wenn diese innerhalb einer formula für arithmetische Berechnungen verwendet werden sollen muss dies innerhalb der I() Funktion erfolgen, siehe R-Hilfe: ?formula.

Mit

lm(y ~ x2 -1)

## 
## Call:
## lm(formula = y ~ x2 - 1)
## 
## Coefficients:
##     x2  
## -2.727

kann eine Regression durch den Ursprung (d.h. ohne Interzept) gerechnet werden.

Als Beispiel legen wir einen data.frame mit dem Namen mydf und sieben Beobachtungen an

mydf <- data.frame(x1 = c(9, 7, 8, 7, 5, 2, 3), 
                   x2 = 1:7, 
                   y  = c(1, 4, 5, 7, 7, 9, 8))

Eine lineare Regression erhalten wir mit

eq1 <- lm(y ~ x1 + x2, data = mydf) 
eq1

## 
## Call:
## lm(formula = y ~ x1 + x2, data = mydf)
## 
## Coefficients:
## (Intercept)           x1           x2  
##      2.0984      -0.0929       1.0757

Für lm Objekte existieren eine Reihe von generischen Funktionen, die sehr nützlich sind (generisch bedeutet sehr vereinfacht, dass sie den Objekttyp erkennen und daraus Parameter übernehmen können)

`print()`	am Bildschirm ausgeben (meist nur innerhalb von Funktionen/Schleifen erforderlich, sonst reicht der Objekt-Name)
`summary()`	kompletten (Regressions-) Output am Bildschirm ausgeben
`coef()`	(oder `coefficients()`) Regressionskoeffizienten extrahieren
`resid()`	(oder `residuals()`) Residuen extrahieren
`fitted()`	gefittete Werte extrahieren
`deviance()`	Quadratsumme der Residuen extrahieren
`vcov()`	Varianz- Kovarianzmatrix der Koeffizienten extrahieren
`plot()`	diagnostische Grafiken für Regression
usw.

Zum Bespiel

summary(eq1)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = mydf)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -1.33802  0.40047  0.41764  1.24902 -0.01249  0.63310 -1.34973 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   2.0984     4.9963   0.420    0.696
## x1           -0.0929     0.4754  -0.195    0.855
## x2            1.0757     0.5742   1.873    0.134
## 
## Residual standard error: 1.215 on 4 degrees of freedom
## Multiple R-squared:  0.8683, Adjusted R-squared:  0.8024 
## F-statistic: 13.19 on 2 and 4 DF,  p-value: 0.01735

Technisch gesehen sind lm Objekte Listen, und mit dem $ Operator kann auf einzelne Elemente der Liste zugegriffen werden, z.B.

eq1$coefficients

## (Intercept)          x1          x2 
##  2.09836066 -0.09289617  1.07572209

Mit der (fast immer nützlichen) str() Funktion erhält man u.a. die Namen der Elemente, auf die zugegriffen werden kann, z.B. str(eq1), oder wenn man nur an deren Namen interessiert ist,

names(eq1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

Verwirrenderweise gibt die auf ein lm Objekt angewandte summary() Funktion ein summary.lm Objekt zurück, in dem auf weitere Elemente zugegriffen werden kann, z.B. auf das Bestimmtheitsmaß

summary(eq1)$r.squared

## [1] 0.8682956

Mit

names(summary(eq1))

##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

erfährt man die Namen dieser Elemente, und mit

summary(eq1)$coefficients

##                Estimate Std. Error    t value  Pr(>|t|)
## (Intercept)  2.09836066  4.9963220  0.4199811 0.6960771
## x1          -0.09289617  0.4753784 -0.1954152 0.8545930
## x2           1.07572209  0.5742413  1.8732927 0.1343186

erhält man ein matrix Objekt mit Koeffizienten, Standardfehler, t-Statistiken und p-Werten. Auf Einzelwerte (oder Vektoren) könnten wir wie üblich bei Matrix Objekten mit den Indizes zugreifen, z.B. auf den p-Wert von x1 (2. Zeile, 4. Spalte) mit der folgenden Funktion zugreifen

summary(eq1)$coefficients[2, 4]

## [1] 0.854593

Dies wird sich später noch öfter als nützlich erweisen.

Der `anova()` Befehl

Mit anova() können mehrere Koeffizienten simultan getestet werden, indem ein kurzes und langes Modell verglichen wird.

Als Beispiel erzeugen wir einen einfachen Datensatz

df <- data.frame(x2 = c(1, 2, 3, 4, 5, 6, 7),
                 x3 = c(8, 9, 7, 3, 4, 1, 0),
                 y =  c(3, 1, 2, 7, 9, 8, 9))
eq_ur <- lm(y ~ x2 + x3, data = df)  ## lang, unrestricted
eq_r  <- lm(y ~ x2, data = df)       ## kurz, restricted

anova(eq_ur, eq_r)

## Analysis of Variance Table
## 
## Model 1: y ~ x2 + x3
## Model 2: y ~ x2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      4 12.070                           
## 2      5 17.393 -1   -5.3224 1.7638 0.2549

In diesem einfachen Beispiel wird nur ein einziger Koeffizient getestet, x3 (da hier nur eine Hypothese getestet wird, könnte dies auch mit einem einfachen t-Test getestet werden, und tatsächlich ist der Wert der F-Statistik gleich dem Quadrat der t-Statistik von x3).

Da der $p$-Wert der F-Statistik (0.2549) größer als 0.05 ist können wir die Nullhypothese, dass der wahre Koeffizient gleich Null ist, nicht verwerfen.

Trickreich wird es mit fehlenden Werten. Nehmen wir an, dass der erste Wert von x2 sowie der dritte und fünfte Wert von x3 fehlen

df$x2[1] <- NA
df$x3[c(3,5)] <- NA
df

##   x2 x3 y
## 1 NA  8 3
## 2  2  9 1
## 3  3 NA 2
## 4  4  3 7
## 5  5 NA 9
## 6  6  1 8
## 7  7  0 9

In diesem Fall beruhen das restringierte und das nicht restringierte Modell auf einer unterschiedlichen Zahl von Beobachtungen, und der anova() Befehl funktioniert nicht mehr!

eq_ur <- lm(y ~ x2 + x3, data = df)  ## lang, unrestricted
eq_r  <- lm(y ~ x2, data = df)       ## kurz, restricted

anova(eq_r, eq_ur)
## Fehler in anova.lmlist(object, ...) : 
##   models were not all fitted to the same size of dataset

Wir könnten mit na.omit(df) alle fehlenden Werte des data.frame entfernen, aber wenn wir einen großen Datensatz mit vielen Variablen haben würden wir vermutlich eine Menge Beobachtungen entfernen, die in nicht benötigten Variablen NA's aufweisen.

Eine nicht sehr elegante, aber funktionierende Lösung ist

eq_ur <- lm(y ~ x2 + x3, data = df)
eq_r  <- lm(y ~ x2, data = na.omit(df[ , all.vars(formula(eq_ur))]))

anova(eq_r, eq_ur)

## Analysis of Variance Table
## 
## Model 1: y ~ x2
## Model 2: y ~ x2 + x3
##   Res.Df    RSS Df Sum of Sq     F  Pr(>F)  
## 1      2 5.1864                             
## 2      1 0.0714  1     5.115 71.61 0.07488 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R Befehle (`functions`)

h.st.

Grundlagen

Eigene Funktionen

Der `lm()` Befehl

Der `anova()` Befehl

R Befehle (functions)

h.st.

Grundlagen

Eigene Funktionen

Der lm() Befehl

Der anova() Befehl

R Befehle (`functions`)

Der `lm()` Befehl

Der `anova()` Befehl