A well-known property of regression models is that they capture the unique contribution of a predictor. By “unique” we mean the effect of the predictor (on the criterion) if the other predictor(s) is/are held constant. A typical classroom example goes along the following lines.

All about storks

1. There’s a correlation between babies and storks. Counties with lots of storks enjoy large number of babies and v.v.

2. However, I have children, I know the storks are not overly involved in that business, so says the teacher (polite laughters in the audience).

3. Once we consider a third variable, call it nature-level, we see that the correlation between storks and babies vanishes.

4. Assume that nature-level and storks are correlated - in nature-rich areas there are a lot of storks. Assume further that nature-level and the amount of babies are correlated - in nature-rich areas people have relatively many babies.

5. That’s why we find a correlation between babies and storks (“spurious” it is called), which is (fully) attributable to the nature-level.

This is all well-known but it helps to see it in real data. That’s the objective of this post.

Tip, sex, and total_bill - two simple bivariate analyses

Let’s have a look at the tips dataset.

data(tips, package = "reshape2")

And load some datajudo packages.

library(tidyverse)

First, simple bivariate assocations for sex-tip and total_bill-tip.

Is there an assocation between sex and tip?

ggplot(tips) +
  ggplot2::aes(x = sex, y = tip) +
  geom_jitter(color = "grey80") +
  geom_jitter(data = tips, aes(color = sex)) +
  stat_summary(fun.y = mean, color = "red", geom = "line", size = 2, group = 1) +
  stat_summary(fun.y = mean, color = "red", geom = "point", size = 5, group = 1, alpha = .5)

tips2 <- select(tips, -sex)

tips %>% 
  ggplot() +
  aes(x = total_bill, y = tip) +
  geom_point(data = tips2, color = "grey80") +
  geom_point(aes(color = sex)) +
  facet_wrap(~sex)

It appears that men do tip slightly more than women do.

tips %>% 
  group_by(sex) %>% 
  summarise(tip_sex_mean = mean(tip),
            tip_sex_md = median(tip))

## # A tibble: 2 x 3
##      sex tip_sex_mean tip_sex_md
##   <fctr>        <dbl>      <dbl>
## 1 Female     2.833448       2.75
## 2   Male     3.089618       3.00

Ok, some 25 Cent difference, in favor of men. Some association is present in the data.

Let’s run a regression on that, just to double-check.

lm(tip ~ sex,
   data = tips)

## 
## Call:
## lm(formula = tip ~ sex, data = tips)
## 
## Coefficients:
## (Intercept)      sexMale  
##      2.8334       0.2562

Ok, confirmed.

Now, what about the association of total_bill and tip?

tips %>% 
  ggplot +
  aes(x = total_bill, y = tip) +
  geom_point()

Oh, clearly, some correlation is present between total_bill and tip.

tips %>% 
  summarise(cor_tip_bill = cor(tip, total_bill))

##   cor_tip_bill
## 1    0.6757341

And let’s have a look at the slope of the regression.

lm(tip ~ total_bill, data = tips)

## 
## Call:
## lm(formula = tip ~ total_bill, data = tips)
## 
## Coefficients:
## (Intercept)   total_bill  
##      0.9203       0.1050

Multivariate analysis

Now let’s do that in one run:

lm(tip ~ total_bill + sex,
   data = tips)

## 
## Call:
## lm(formula = tip ~ total_bill + sex, data = tips)
## 
## Coefficients:
## (Intercept)   total_bill      sexMale  
##     0.93328      0.10523     -0.02661

Notice that the slope of the predictor total_bill has not changed much. But sex did!

So, the unique contribution of sex is very small, about 2 Cent, and in favor of women when considered in junction with total_bill!

Let’s visualize the multivariate appraoch

First, for easier visualization, let’s bin total_bill in say, 5, bins.

tips$bill_bins <- cut_interval(tips$total_bill, n = 5)

Now let’s plot all three variables in one go.

tips %>% 
  ggplot +
  aes(x = bill_bins, 
      y = tip,
      fill = sex) +
  geom_boxplot() +
  geom_point(alpha = .5)

We see that females do tip more than men when the 5 bins of total_bill are considered each in isolation. For 4 of the 5 bins we find that the median of the women is higher than that of the men. Only in the 5th, ie., highest bin, we found that men tip higher. However, this bin is slightly populated, so we may negligate it. In sum, looking at each bin of total_bill - women tip more than men do!

Depicted differently:

tips %>% 
  ggplot() +
  aes(x = total_bill, y = tip) +
  geom_point(data = tips2, color = "grey80") +
  geom_point(aes(color = sex)) +
  facet_wrap(~sex)

Dichotomize total_bill

Maybe even more vivid is we we do not bin total_bill in 5 bins, but just in two.

tips$bill_bins <- cut_interval(tips$total_bill, n = 2)

Again, let’s plot the three variables.

tips %>% 
  ggplot +
  aes(x = bill_bins, 
      y = tip,
      fill = sex) +
  geom_boxplot() +
  geom_point(alpha = .5)

We find that women seem to tip slightly more than men do; at least not less.

A Simpson case

What happens can be considered a Simpson’s paradox, based on the fact that men spend more than women do (higher total_bill). Even though they tip less or equally to women it appears they tip “in sum” (not considering total_bill) more simply because total_bill and tip are related.

Let’s dichotomize total_bill and see which sex group (females or males) occur more frequently.

tips %>% 
  mutate(bill_di = ntile(total_bill, n = 2),
         bill_di_f = if_else(bill_di == 1, "low", "hi")) %>% 
  count(sex, bill_di_f) %>% 
 # gather(bill_group, n, -n) %>% 
  ggplot() +
  aes(x = sex, fill = bill_di_f, y = n) +
  geom_col(position = "fill")

We see that the in the male group, high values occurr more frequently.

Path model

sex_male has a positive effect on total_bill; and total_bill in turn on tip. But sex_male does not exert much of an influence on tip anymore, after controling for total_bill (in fact, the positive effect turned to a weak negative impact).

library(DiagrammeR)
mermaid("
graph LR
  total_bill-->tip
  sex_male-->total_bill
  sex_male---tip
 ")

Conclusion

That kind of “swapping” of the regression coefficient is a frequent observation in (multiple) regression. We may conclude that the multivariate analyses with 3 variables is more adequate than the bivariate analyses with one predictor only.

However, a caveat is that we can never be sure that this swapping process is now complete. Maybe we are missing out some more variable that will again make the coefficients swap. Alas, this we cannot know in observational reseach. That’s why observational data should not be overinterpreted.

Of course, even if commonly speak of “effects” this language should never be taken as indicative of causal effects. Causal effects is nothing that can firmly be established by statistics. Causality claims rest of more logical and substantive argument.

Nicht jeder liebt Datenanalyse und Statistik… in gleichem Maße. Das ist zumindest meine Erfahrung aus dem Unterricht . Crashkurse zu R sind vergleichbar zu Crahskursen zu Französisch - kann man machen, aber es sollte die Maxime gelten “If everything else fails”.

Dieser Crashkurs ist für Studierende oder Anfänger der Datenanalyse gedacht, die in kurzer Zeit einen verzweifelten Versuch … äh … einen grundständigen Überblick über die Datenanalyse erwerben wollen.

Ok, also los geht’s im Schweingalopp durch die Datenanalyse .

Software

Bevor wir uns die Schritte näher anschauen, ein paar Worte zur Software.

Programme

Wir brauchen zwei Programme.

• R
• RStudio

Wenn R installiert ist, dann findet RStudio R auch direkt. Wenn alles läuft, sieht es etwa so aus:

Erweiterungen

Außerdem brauchen wir diese Pakete für R:

• tidyverse
• mosaic

Um einen Befehl zu verwenden, der nicht im Standard-R, sondern in einer Erweiterung von R (“Paket”) wohnt, müssen sie dieses Paket erst starten (laden). Dazu können Sie den Befehl library verwenden. Oder Sie klicken den Namen des Pakets hier an:

Um ein Paket zu laden, muss es installiert sein. Klicken Sie auf den Button “Install” unter dem Reiter “Packages” in RStudio:

Starten Sie jetzt die beiden Erweiterungen (per Klick oder mit folgenden Befehlen).

library(mosaic)
library(tidyverse)

Über sieben Brücken musst Du gehen

Lizenz: André D Conrad, CC BY SA 3.0 De, https://de.wikipedia.org/wiki/Peter_Maffay#/media/File:Peter_Maffay.jpg

Man kann (wenn man will) die Datenanalyse in sieben fünf Brücken oder Schritte einteilen, angelehnt dem Song von Peter Maffay “Über sieben Brücken musst du gehen”.

Brücke 1: Daten einlesen

Der einfachste Weg, Daten einzulesen, ist über den Button “Import Dataset” in RStudio.

So lassen sich verschiedene Formate - wie XLS(X) oder CSV - importieren.

Beim Importieren von CSV-Dateien ist zu beachten, dass R davon von us-amerikanisch formatierten CSV-Dateien ausgeht. Was heißt das? Das bedeutet, das Spaltentrennzeichen (delimiter) ist ein Kommma ,. Deutsch formatierte CSV-Dateien, wie sie ein deutsch-eingestellltes Excel ausgibt, nutzen aber ein Semikolon ; (Strichpunkt) als Spaltentrennzeichen.

Ggf. müssen Sie also in der Import-Maske von RStudio als delimiter ein semicolon auswählen.

Alternativ können Sie natürlich eine XLS- oder XLSX-Datei importieren.

Am einfachsten ist es, XLSX-Dateien zu importieren.

tidy data - Tabellen in Normalform

Damit Sie in R vernünftig mit Ihren Daten arbeiten können, sollten die Daten “tidy” sein, d.h. in Normalform. Was ist Normalform? Betrachten Sie folgende Abbildung:

Die goldene Regel der Normalform einer Tabelle lautet also:

In jeder Zeile steht eine Beobachtung (z.B. Person). In jeder Spalte eine Variable (z.B. Geschlecht). In der ersten Zeile stehen die Spaltennamen, danach folgen die Werte. Sonst steht nichts in der Tabelle.

Falls Ihre Daten nicht in Normalform sind, sollten Sie diese zunächst in Normalform bringen.

Der einfachste Weg (von der Lernkurve her betrachtet, nicht vom Zeitaufwand), Daten in Normalform zu bringen, ist sie in Excel passend umzubauen.

Beispiel für Daten in Nicht-Normalform

Sie denken, dass Ihre Daten immer/auf jeden Fall in Normalform sind? Dann scheuen Sie sich mal dieses Bild an:

Textkodierung in UTF-8

Falls Sie RStudio oder ein beliebiger Texteditor irgendwann fragt, wie die Textdatei kodiert sein soll, wählen Sie immer “UTF-8”. Googeln Sie bei Interesse, was da bedeutet.

Schritt 2: Aufbereiten

Der Schritt des Aufbereitens ist häufig der zeitintensivste Schritt. In diesem Schritt erledigen Sie alles, bevor Sie zu den “coolen” oder fortgeschrittenen Analysen kommen. Z.B.

• prüfen auf Fehler beim Daten einlesen (und korrigieren)
• Spaltennamen korrigieren
• Daten umkodieren
• Fehlende Werte verarzten
• Komische Werte prüfen
• Daten zusammenfassen
• Zeilenmittelwerte bilden
• Logische Variablen bilden

Auf Fehler prüfen

Ein häufiger Fehler ist, dass die Daten nicht richtig eingelesen werden. Zum Beispiel werden die Spaltentrennzeichen nicht richtig erkannt. Das kann dann so aussehen.

Unter “delimiter” in der Maske können Sie das Trennzeichen anpassen.

“Deutsche” CSV-Dateien verwenden als Dezimaltrennzeichen ein Komma; englisch-formatierte CSV-Dateien hingegen einen Punk. R geht per Default von englisch-formatierten CSV-Dateien aus. Importieren Sie eine deutsch-formatierte CSV-Datei, müssen Sie das Dezimaltrennzeichen von Hand ändern; es wird nicht automatisch erkannt.

Unter “locale” können Sie das Dezimaltrennzeichen anpassen.

Spaltennamen korrigieren

Spaltennamen müssen auch “tidy” sein. Das heißt in diesem Fall:

• keine Leerzeichen
• keine Sonderzeichen (#,ß,ä,…)
• nicht zu lang, aber trotzdem informativ

Spaltennamen sollten nur Buchstaben (ohne Umlaute) und Ziffern enthalten.

Für Textdaten in den Spalten sind diese Regeln auch sinnvoll.

Umkodieren

Gerade bei der Analyse von Fragebogendaten ist es immer wieder nötig, Daten umzukodieren. Klassisches Besispiel: Ein Item ist negativ kodiert. Zum Beispiel das Item “Ich bin ein Couch-Potato” in einem Fragebogen für Extraversion.

Nehmen wir an, das Item “i04” hat die Werte 1 (“stimme überhaupt nicht zu”) bis 4 (“stimme voll und ganz zu”). Kreuzt jemand das Couch-Potato-Item mit 4 an, so sollte er nicht die maximale Extraversion-Punktzahl (4), sondern die minimale Extraversion-Punktzahl (1) erhalten. Also

1 --> 4 2 --> 3 3 --> 2 4 --> 1

Am einfachsten ist dies zu bewerkstelligen mit folgendem R-Befehl:

meine_tabelle$i04_r <- 5 - meine_Tabelle$i04

Rechnet man 5-i04 so kommt der richtige, “neue” Wert heraus.

Zur Erinnerung:

• $ ist das Trennzeichen zwischen Tabellennamen und Spaltenname.
• <- ist der Zuweisungsbefehl. Wir definieren eine neue Spalte mit dem Namen i04_r. Das r soll stehen für “rekodiert”, damit wir wissen, dass in dieser Spalte die umkodierten Werte stehen.

Fehlende Werte

Der einfachste Umgang mit fehlenden Werten ist: nichts machen. Denken Sie nur daran, dass viele R-Befehle von Natur aus nervös sind - beim Anblick von fehlenden Werten werden sie panisch und machen nix mehr. Zum Beispiel der Befehl mean. Haben sie fehlende Werte in ihren Daten, so verwenden Sie den Parameter na.rm = TRUE. na steht für “not available”, also fehlende Werte. rm steht für “remove”. Also mean(meine_tabelle$i04_r, na.rm = TRUE).

Der R-Befehl summary zeigt Ihnen an, ob es fehlende Werte gibt:

summary(meine_daten).

Komische Werte

Hat ein Spaßvogel beim Alter 999 oder -1 angegeben, kann das Ihre Daten ganz schön verhageln. Prüfen Sie die Daten auf komische Werte. Der einfachste Weg ist, sich die Daten in Excel anzuschauen. Cleverer ist noch, sich Zusammenfassungen auszugeben, wie der kleinste oder der größte Wert, oder der Mittelwert etc., und dann zu schauen, ob einem etwas spanisch vorkommt. Diagramme sind ebenfalls hilfreich.

Logische Variablen bilden

Sagen wir, uns interessiert welches Auto mehr als 200 PS hat; wir wollen Autos mit mehr als 200 PS vergleichen (“Spass”) mit schwach motorisierten Autos (“Kruecke”). Wie können wir das (einfach) in R erreichen? Logische Variablen sind ein einfacher Weg.

mtcars$Spass <- mtcars$hp > 200

Dieser Befehl hat eine Spalte (Variable) in der Tabelle mtcars erzeugt, in der TRUE steht, wenn das Auto der jeweiligen Spalte die Bedingung (hp > 200) erfüllt. Schauen Sie nach.

favstats(mtcars$Spass)

## Warning in fav_stats(x, ..., na.rm = na.rm): Auto-converting logical to
## numeric.

##  min Q1 median Q3 max    mean        sd  n missing
##    0  0      0  0   1 0.21875 0.4200134 32       0

Daten zusammenfassen: Deskriptivstatistik

Deskriptive Statistik ist letztlich nichts anderes, als Daten geschickt zusammenzufassen.

Praktisch sind es meistens Spalten, die zu einen Zahl zusammengefasst werden.

Schauen wir uns das mal mit echten Daten an. Der Datensatz “mtcars” ist schon in R eingebaut, so dass wir in nicht extra laden müssen. Ganz praktisch.

summary(mtcars)

##       mpg             cyl             disp             hp       
##  Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0  
##  1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5  
##  Median :19.20   Median :6.000   Median :196.3   Median :123.0  
##  Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7  
##  3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0  
##  Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0  
##       drat             wt             qsec             vs        
##  Min.   :2.760   Min.   :1.513   Min.   :14.50   Min.   :0.0000  
##  1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000  
##  Median :3.695   Median :3.325   Median :17.71   Median :0.0000  
##  Mean   :3.597   Mean   :3.217   Mean   :17.85   Mean   :0.4375  
##  3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000  
##  Max.   :4.930   Max.   :5.424   Max.   :22.90   Max.   :1.0000  
##        am              gear            carb         Spass         cyl_f 
##  Min.   :0.0000   Min.   :3.000   Min.   :1.000   Mode :logical   4:11  
##  1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000   FALSE:25        6: 7  
##  Median :0.0000   Median :4.000   Median :2.000   TRUE :7         8:14  
##  Mean   :0.4062   Mean   :3.688   Mean   :2.812   NA's :0               
##  3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000                         
##  Max.   :1.0000   Max.   :5.000   Max.   :8.000

Hilfe zu diesem Datensatz bekommen Sie so:

help(mtcars)

mit help(Befehl) bekommt man Hilfe zu einem Befehl oder einem sonstigen Objekt (z.B. Datensatz).

Numerische Variablen

Der einfachste Weg, um Deskriptivstatistik für eine numerische Variable auf einen Abwasch zu erledigen ist dieser Befehl:

favstats(mtcars$mpg)

##   min     Q1 median   Q3  max     mean       sd  n missing
##  10.4 15.425   19.2 22.8 33.9 20.09062 6.026948 32       0

Der Befehl favstats lässt auch Subgruppenanalysen zu, z.B. um Männer und Frauen zu vergleichen:

favstats(mpg ~ cyl, data = mtcars)

##   cyl  min    Q1 median    Q3  max     mean       sd  n missing
## 1   4 21.4 22.80   26.0 30.40 33.9 26.66364 4.509828 11       0
## 2   6 17.8 18.65   19.7 21.00 21.4 19.74286 1.453567  7       0
## 3   8 10.4 14.40   15.2 16.25 19.2 15.10000 2.560048 14       0

Dabei ist mpg die Variable, die sie vergleichen wollen (Spritverbrauch); cyl die Gruppierungsvariable (Anzahl der Zylinder). Gruppierungsvariable bedeutet hier, dass den Spritverbrauch zwischen 4,6 und 8-Zylindern vergleichen wollen.

Nominale Variablen

Eine Häufigkeitstabelle für eine nicht-metrische Variable lässt über den Befehl table erstellen.

Aber zuerst wandeln wir eine unserer metrischen Variablen in eine nominale um:

my_cars <- mtcars

my_cars$cyl_f <- factor(mtcars$cyl)
my_cars$Spass <- my_cars$hp > 200

factor wandelt eine metrische (numerische) Variable in einen nominale Variable um; nominale Variablen nennt man in R auch “Faktor”. Um den “originalen” Datensatz nicht zu überschreiben, packen wir alles in eine Kopie mit dem Namen my_cars.

table(my_cars$cyl_f)

## 
##  4  6  8 
## 11  7 14

Aha, 14 Autos in unserer Tabelle haben also 8 Zylinder.

Brechen wir jetzt die Auszählunng noch mal auf in “Spass” vs. “Krücke”.

table(my_cars$cyl_f, my_cars$Spass)

##    
##     FALSE TRUE
##   4    11    0
##   6     7    0
##   8     7    7

Zeilenmittelwerte bilden

Bei Umfragen kommt es häufig vor, dass man Zeilenmittelwerte bildet. Wieso? Man möchte z.B. in einer Mitarbeiterbefragung den “Engagementwert” jedes Beschäftigten wissen (klingt einfach gut). Dazu addiert man die Werte jedes passenden Items auf. Diese Summe teilen Sie durch die Anzahl der Spalten

Zeilenmittelwerte bilden Sie am einfachsten in Excel.

In R können Sie Zeilen einfach mit dem + Zeichen addieren:

meine_tabelle$zeilenmittelwert <- (meine_tabelle$item1 + meine Tabelle$item2) / 2

Schritt 3: Visualisieren

Ein Bild sagt bekanntlich mehr als 1000 Worte. Betrachten Sie dazu “Anscombes Quartett”:

Diese vier Datensätze sehen ganz unterschiedlich aus, nicht wahr? Aber ihre zentralen deskriptiven Statistiken sind praktisch gleich! Ohne Diagramm wäre uns diese Unterschiedlichkeit nicht (so leicht) aufgefallen!

Zur Visualisierung empfehle ich das R-Paket ggplot2. Es wird mit geladen, wenn Sie tidyverse laden.

Es gib einen recht einfachen Befehl in diesem Paket, qplot. Bei qplot geben Sie folgende Parameter an:

• X-Achse (x): Welche Variable (Spalte in ihrer Tabelle) soll auf der X-Achse stehen?
• Y-Achse (y): Welche Variable (Spalte in ihrer Tabelle) soll auf der Y-Achse stehen?
• Daten (data): Wie heißt ihre Datentabelle?
• Geom (geom): Welche Art von Bildchen wollen Sie malen (z.B. Punkte, Boxplots, Linien, Histogramm…)`

Darüber hinaus verkrafter der Befehl noch viele andere Schnörkel, die wir uns hier sparen. Interessierte können googeln… Es ist ein sehr mächtiger Befehl, der sehr ansprechende Diagramme erzeugen kann.

Probieren wir’s!

qplot(x = hp,
      y = mpg,
      geom = "point",
      data = mtcars)

Easy, oder?

Ein anderes Geom:

qplot(x = factor(cyl),
      y = mpg,
      data = mtcars,
      geom = "boxplot")

Beachten Sie, dass qplot nur dann mehrere Boxplots zeichnet, wenn auf der X-Achse eine nominal skalierte Variable steht.

Eine metrische Variable wandeln Sie in eine nominal skalierte Variable um mit dem Befehl factor(meine_metrische_variable).

Oder mal nur eine Variable:

qplot(x = hp,
      data = mtcars,
      geom = "histogram")

Geben wir keine Y-Variable an, nimmt qplot eigenständig die Häufigkeit pro X-Wert!

Tauschen Sie mal “histogram” mit “density”!

Schritt 4: Modellieren

Modellieren hört sich kompliziert an. Für uns hier heißt es vor allem ein inferenzstatistisches Verfahren anzuwenden.

Das Schweizer Taschenmesser und den Modellierungsverfahren ist die Regressionsanalyse. Man kann sie für viele Zwecke einsetzen.

Wann welchen Test?

Es gibt in vielen Lehrbüchern Übersichten zur Frage, wann man welchen Test rechnen soll. Googeln hilft hier auch weiter. Eine Übersicht findet man hier: http://www.methodenberatung.uzh.ch/de/datenanalyse.html.

Wie heißt der jeweilige R-Befehl?

Wenn man diese Befehle nicht häufig verwendet, ist es schwierig, sie auswändig zu wissen. Googeln Sie. Eine gute Übersicht findet sich hier: http://r-statistics.co/Statistical-Tests-in-R.html.

Regression

Weil die Regression so praktisch ist, hier ein Beispiel.

lm(mpg ~ cyl, data = mtcars)

## 
## Call:
## lm(formula = mpg ~ cyl, data = mtcars)
## 
## Coefficients:
## (Intercept)          cyl  
##      37.885       -2.876

lm heißt “lineares Modell” - weil man bei der (normalen) Regression eine Gerade in die Punktewolke der Daten legt, um den Trend zu abzuschätzen. Als nächstes gibt man die “Ziel-Variable” (Output) an, hier mpg. Dann kommt ein Kringel ~ gefolgt von einer (mehr) Input-Variablen (Prädiktoren, UVs). Schließnoch muss noch die Datentabelle erwähnt werden.

Das Ergbebnis sagt uns, dass pro Zylinder die Variable mpg um knapp 3 Punkte sinkt. Also: Hat die Karre einen Zylinder mehr, so kann man pro Galone Sprit 3 Meilen weniger fahren. Immer im Schnitt, versteht sich. (Und wenn die Voraussetzungen erfüllt sind, aber darum kümmern wir uns jetzt nicht.)

Allgemein:

lm(output ~ input, data = meine_daten)

Easy, oder?

Man kann auch mehrere Prädiktoren anführen:

lm(mpg ~ cyl + hp, data = mtcars)

## 
## Call:
## lm(formula = mpg ~ cyl + hp, data = mtcars)
## 
## Coefficients:
## (Intercept)          cyl           hp  
##    36.90833     -2.26469     -0.01912

Dazu werden die durch + getrennt. Die Ergebnisse zeigen uns, dass die PS-Zahl (´hp´) kaum Einfluss auf den Spritverbrauch haut. Genauer: Kaum zusätzlichen Einfluss auf den Spritverbrauch hat. Also Einfluss, der über den Einfluss hinausgeht, der schon durch die Anzahl der Zylinder erklärt werden würde. Es ist also praktisch wurscht, wie viel PS das Auto hat, wenn man den Verbrauch schätzen will - Hauptsache, man weiß die Anzahl der Zylinder.

Vorhersagen

Man kann die Regression nutzen, um Vorhersagen zu treffen. Sagen wir, unser neuer Lamborghini hat 400 PS und 12 Zylinder. Wie groß ist wohl der Spritverbrauch laut unserem Regressionsmodell?

Als Vorbereitung speichern wir unser Regressionsmodell in einer eigenen Variablen:

mein_lm <- lm(mpg ~ cyl + hp, data = mtcars)

Dazu nimmt man am besten den Befehl predict, weil wir wollen eine Vorhersage treffen:

predict(mein_lm, data.frame(cyl = 12, hp = 400))

##        1 
## 2.083329

Aha. Pro Gallone kämmen wir 2 Meilen. Schluckspecht.

Visualisierung des Modells

Schauen wir uns das Modell mal an, damit es nicht so theoretisch ist. Wie gesagt, ein Regressionsmodell ist nichts anderes als eine Trendgerade in einer Punktewolke.

Gerade sind durch Achsenabschnitt (engl. “intercept”) f(x=0) und Steigung (englisch: “slope”) definiert. Der Steigungswert ist die Zahl, die der lm-Befehl für den Prädiktor ausspuckt

qplot(x = cyl,
      y = mpg,
      data = mtcars,
      geom = "point") +
  geom_abline(slope = -2.9, intercept = 38)

Schritt 5: Kommunizieren

Kommunizieren soll sagen, dass Sie Ihre Ergebnisse anderen mitteilen - als Student heißt das häufig in Form einer Seminararbeit an den Dozenten.

Einige Hinweise:

• Geben Sie nicht alle Ergebnisse heraus. Ihre Fehler müssen niemanden interessieren.
• Die wesentlichen Ergebnisse kommen in den Hauptteil der Arbeit. Interessante Details in den Anhang.
• Der Mensch ist ein Augentier. Ihr Gutachter auch. Achten Sie auf optisch ansprechende Darstellung; schöne Diagramme helfen.
• Dozenten achten gerne auf formale Korrektheit. Das Gute ist, dass dies relativ einfach sicherzustellen ist, da auf starren Regeln basierend.

Tabellen und Diagramme

Daten kommunizieren heißt praktisch zumeist, Tabellen oder Diagramme zu erstellen. Meist gibt es dazu Richtlinien von Seiten irgendeiner (selbsterannten) Autorität wie Dozenten oder Fachgesellschaften. Zum Beispiel hat die APA ein umfangreiches Manual zum Thema Manuskripterstellung publiziert; die deutsche Fachgesellschaft der Psychologie entsprechend. Googeln Sie mal, wie in ihren Richtlinien Tabellen und Diagramme zu erstellen sind (oder fragen Sie Ihren Gutachter).

Für Fortgeschrittene: RMarkdown

Wäre das nicht cool: Jegliches Formatieren wird automatisch übernommen und sogar so, das ses schick aussieht? Außerdem wird Ihr R-Code und dessen Ergebnisse (Tabellen und Diagramme oder reine Zahlen) automatisch in Ihr Dokument übernommen. Keine Copy-Paste-Fehler mehr. Keine händisches Aktualisieren, weil Sie Daten oder die vorhergehende Analyse geändert haben. Hört sich gut an? Probieren Sie mal RMarkdown aus :one_up:.

Fazit

Glück auf

Hier finden Sie einen Überblick an einige “Lieblings-R-Befehle”.

One way to dig into some topic such as data analysis is just-doing, trial and error. Another way is reading blogs; a fruitful avenue in my experience. However, the classical way of reading some good book is all but outdated.

Here are some recommendations of books I found helpful as a starter (books in English and German).

R for Data Science

Grolemund, G., & Wickham, H. (2016). R for Data Science. O’Reilly Media, Incorporated. Retrieved from https://books.google.de/books?id=aZRYrgEACAAJ

Hadley Wickham did it again. This book is a great ressource for building up a solid base in data analysis using R. What is special is that, of course, the philosophy of “tidy data” is a pivotal path throughout the book. That comes as no surprise as Hadleyi is the major prophet of this way of thinking (a great way by the way). Second, the book is “modern” in the sense as newly developed, but thoughtfully devised and already quite popular approaches are put forward. Note that the book is not intended to go into all the depth of the R languages, neither programming, neither statistical learning. Most notably whereas the title indicates a salty dose of statistics, this association does not come true. Statistical modeling is a rather small topic in the book.

Introduction to Statistical Learning

James, G., Witten, D., Hastie, T., & Tibishirani, R. (2013). An Introduction to Statistical Learning (Vol. 112). New York City: Springer. DOI: http://doi.org/10.1007/978-1-4614-7138-7

A classical introduction, aimed at not-so mathematically inclined (or trained) learners. Important topics on the principles of statistical learnings are covered as well as widely used algorithms. R serves as the computational environment for all models/ computations. Datastes, exercised, and case studies (as well as an R package) comes bundled with the book. A great start into a deep data dive. Great extra: The PDF of the book is freely available in the web.

R - Einführung durch angewandte Statistik

Hatzinger, R., Hornik, K., & Nagel, H. (2011). R: Einführung durch angewandte Statistik. Hallbergmoos: Pearson Studium. Retrieved from http://books.google.de/books?id=jmCuEY07s7UC

This book is intended for (native) German speaking audiences, especially in college courses and the like. The book offers a broad perspective on “typical” statistics using R. The curriculum of the book very well reflects typical college curricula; a notion that may find widespread inclinement throughout teachers and professors. However, the book is not modern in the sense that more recent R approaches such as the “tidyverse” are discussed in the book. Modern critical discussion on the p-value, presentation of alternatives, reflection on open science, and more on statistical learning would have helped to render a good book great, in my opinion. However, to get a grip on standard statistics using R, this book is a valuable resource.

knitr::opts_chunk$set(fig.align = "center",
                      out.width = "70%",
                      fig.asp = .61)

Every now and then, random numbers come in handy to demonstrate some statistical behavior. Of course, well-known appraoches are rnorm and friends. These functions are what is called pseudo random number generators, because they are not random at all, strictly speaking, but determined by some algorithm. An algorithm is a sort of creature that is 100% predictable once you know the input (and the details of the algorithm). Pseudo random numbers are useful ‘cause you know what you will get; reproducibility ensured.

However, at times it is useful to get make use of true random numbers, coming from atmospheric noise, for example. That’s were the R package random by Dirk Eddelbuettel comes into play.

Technical details are beyond the scope of this post; here we just grap some random numbers to simulate or demonstrate some research setting.

For that end, suppose we conducted an experiment with three groups. Our hypothesis was that stories plus images of persons in despair will induce more pro social behavior (ie. amount of donation) compared to presenting stats about suffering people.

Say we came up with three experimental groups:

1. Statistics only (“100,000 are at the brink of famine”)
2. Story plus image of child (“Rodia from Sudan has not eaten since a week”. Picture added.)
3. Combination of 1. and 2.

The basis of this experiment is from this paper:

Slovic, P. (2007). If I look at the mass I will never act: Psychic numbing and genocide. Judgment and Decision Making, 2(2), 79–95. https://doi.org/10.1007/978-90-481-8647-1

For a presentation, I wanted to demonstrate simulated data.

That’s what I did:

First, load some packages.

library(tidyverse)
library(random)

Next, simulate true random data.

Note that the data is uniformally distributed, more or less.

raw <- random::randomNumbers(n = 150, min = 1, max = 10, col = 3) %>% 
  as_tibble

Then, we tidy up the dataframe:

df <- raw %>% 
  mutate(V2 = V2 + 5,
         V3 = V3 + 10) %>% 
  gather %>% 
  mutate(group = recode(key, 
                        V1 = "stat", 
                        V2 = "stat+img",
                        V3 = "img"))

… and plot it:

ggplot(df) +
  aes(x = group, y = value, color = group) +
  geom_boxplot() +
  geom_jitter() +
  theme_minimal() +
  scale_color_brewer(palette = "Set1")

In order sort the groups to our desire, we can make use of factors.

df$group <- factor(df$group, levels = c("stat", "stat+img", "img"))

And plot again, now with groups sorted:

ggplot(df) +
  aes(x = group, y = value, color = group) +
  geom_boxplot() +
  geom_jitter() +
  theme_minimal() +
  scale_color_brewer(palette = "Set1")

The logistic regression is an incredible useful tool, partly because binary outcomes are so frequent in live (“she loves me - she doesn’t love me”). In parts because we can make use of well-known “normal” regression instruments.

But the formula of logistic regression appears opaque to many (beginners or those with not so much math background).

Let’s try to shed some light on the formula by discussing some accessible explanation on how to derive the formula.

Plotting the logistics curve

Let’s have a look at the curve of the logistic regression.

For example, let’s look at the question of whether having extramarital affairs is a function of marital satisfaction.

library(tidyverse)
data(Affairs, package = "AER")

Affairs$has_affairs <- if_else(condition = Affairs$affairs > 0, 1, 0)

Affairs %>% 
  ggplot() +
  aes(x = rating,y = has_affairs) +
  geom_jitter() +
  geom_smooth(aes(y = has_affairs), method = "glm", method.args = list(family = "binomial"), se = FALSE)

Hm, the curve does not pop out to vividly. Let’s have a look at some other data, why not the Titanic disaster data. Can the survival chance plausible put as a function of the cabin fare?

data(titanic_train, package = "titanic")

titanic_train %>% 
  ggplot() +
  aes(x = Fare,y = Survived) +
  geom_point() +
  geom_smooth(aes(y = Survived), method = "glm", method.args = list(family = "binomial"), se = FALSE)

Hm, maybe better to look at the curve in general.

logist <- function(x){
  y = exp(x) / (1 + exp(x))
}

p1 <- ggplot(data_frame())

p1 + stat_function(aes(-5:5), fun = logist) + xlab("x")

Ok, better.

Functional form

It is well-known that the fucntional form of the logictic regression curve is

where \(e\) is Euler’s number (2.718…) and \(t\) can be any linear combination of predictors such as \(b0 + b1x\). \(Y=1\) indicates that the event in question has occured (eg., “survived”, “has_affairs”).

Assume that \(t\) is \(b0 + b1x\) then

Now what? Well, we would to end up with the “typical” formula of the logistic regression, something like:

where \(L\) is the Logit, i.e.,

Deriving the formula

Ok, in a first step, let’s take our \(p(Y=1) = f(t)\) and divide by the probability of the complementary event. If the probability of event \(A\) is \(p\), the the probability of \(not-A\) is \(1-p\). Thus

So wat did we do? We have just replaced \(f(t)\) by \(\frac{e^t}{1e^t}\), and have therey computed the odds.

Next, we multiply the equation by \(\frac{1+e^t}{1+e^t}\) (which is the neutral element, 1), yielding.

In other words, the denominator of the numerator “wandered” down to the denominator.

Now, we can simplify the denominator a bit:

Simplifying the denominator further

But the denominator simplifies to $1$, as can be seen here

so the final solution is

\(=e^t\).

Ok, great, but what does this solution tells us? It tells us the that the odds simplify to \(e^t\).

Now, let’s take the (natural) logarithm of this expression.

by the rules of exponents algebra.

But \(t = b0 + b1x\).

In sum

The left part of the previous equation is called the logit which is “odds plus logarithm” of \(f(t)\), or rather, more precisely, the logarithm of the odss of \(p/(1-p)\).

Looking back, what have we gained? We now know that if we take the logit of any linear combination, we will get the logistic regression formula. In simple words: “Take the normal regression equation, apply the logit \(L\), and you’ll get out the logistic regression” (provided the criterion is binary).

\(L(t) = ln\left( \frac{f(t)}{1-f(t)}\right) = b0 + b1x\).

The formula of the logistic regression is similar in the “normal” regression. The only difference is that the logit function has been applied to the “normal” regression formula.

The linearity of the logit helps us to apply our standard regression vocabulary: “If X is increased by 1 unit, the logit of Y changes by b1”. Just insert “the logit”; the rest of the sentence is the normal regression parlance.

Note that the slope of the curve is not linear, hence b1 is not equal for all values of X.