Hilfe! Mein R startet nicht! Mein R startet zwar, tut aber nicht so, wie ich will. Sicherlich hat es sich (wieder einmal) gegen mich verschworen. Wahrscheinlich hilft nur noch Verschrotten… Bevor Sie zum äußersten schreiten, hier einige Tipps, die sich bewährt haben.

Lösungen, wenn R nicht (richtig) läuft

  • AEG: Aus. Ein. Gut. Starten Sie den Rechner neu. Gerade nach Installation neuer Software zu empfehlen.

  • Sehen Sie eine Fehlermeldung, die von einem fehlenden Paket spricht (z.B. “Package ‘Rcpp’ not available”) oder davon spricht, dass ein Paket nicht installiert werden konnte (z.B. “Package ‘Rcpp’ could not be installed” oder “es gibt kein Paket namens ‘Rcpp’” oder “unable to move temporary installation XXX to YYY”), dann tun Sie folgendes:

    • Schließen Sie R und starten Sie es neu.
    • Installieren Sie das oder die angesprochenen Pakete mit install.packages("name_des_pakets", dependencies = TRUE) oder mit dem entsprechenden Klick in RStudio.
    • Starten Sie das entsprechende Paket mit library(paket_name).
  • Gerade bei Windows 10 scheinen die Schreibrechte für R (und damit RStudio oder RComannder) eingeschränkt zu sein. Ohne Schreibrechte kann R aber nicht die Pakete (“packages”) installieren, die Sie für bestimmte R-Funktionen benötigen. Daher schließen Sie R bzw. RStudio und suchen Sie das Icon von R oder wenn Sie RStudio verwenden von RStudio. Rechtsklicken Sie das Icon und wählen Sie “als Administrator ausführen”. Damit geben Sie dem Programm Schreibrechte. Jetzt können Sie etwaige fehlende Pakete installieren.

  • Ein weiterer Grund, warum R bzw. RStudio die Schreibrechte verwehrt werden könnnten (und damit die Installation von Paketen), ist ein Virenscanner. Der Virenscanner sagt, nicht ganz zu Unrecht, “Moment, einfach hier Software zu installieren, das geht nicht, zu gefährlich”. Grundsätzlich gut, in diesem Fall unnötig. Schließen Sie R/RStudio und schalten Sie dann den Virenscanner komplett (!) aus. Öffnen Sie dann R/RStudio wieder und versuchen Sie fehlende Pakete zu installieren.

  • Läuft der RCommander unter Mac nicht, dann prüfen Sie, ob Sie X11 installiert haben. X11 muss installiert sein, damit der RCommander unter Mac läuft.

  • Die “app nap” Funktion beim Mac kann den RCommander empfindlich ausbremsen. Schalten Sie diese Funktion aus z.B. im RCommander über Tools - Manage Mac OS X app nap for R.app.

Allgemeine Hinweise zur Denk- und Gefühlswelt von R

  • Wenn Sie RStudio starten, startet R automatisch auch. Starten Sie daher, wenn Sie RStudio gestartet haben, nicht noch extra R. Damit hätten Sie sonst zwei Instanzen von R laufen, was zu Verwirrungen (bei R und beim Nutzer) führen kann.

  • Ein neues R-Skript im RStudio können Sie z.B. öffnen mit File-New File-R Script.

  • R-Skripte können Sie speichern (File-Save) und öffnen.

  • R-Skripte sind einfache Textdateien, die jeder Texteditor verarbeiten kann. Nur statt der Endung txt, sind R-Skripte stolzer Träger der Endung R. Es bleibt aber eine schnöde Textdatei.

  • Bei der Installation von Paketen mit install.packages("name_des_pakets") sollte stets der Parameter dependencies = TRUE angefügt werden. Also install.packages("name_des_pakets", dependencies = TRUE). Hintergrund ist: Falls das zu installierende Paket seinerseits Pakete benötigt, die noch nicht installiert sind (gut möglich), dann werden diese sog. “dependencies” gleich mitinstalliert (wenn Sie dependencies = TRUE setzen).

  • Hier finden Sie weitere Hinweise zur Installation des RCommanders: http://socserv.socsci.mcmaster.ca/jfox/Misc/Rcmdr/installation-notes.html.

  • Sie müssen online sein, um Packages zu installieren.

  • Die “app nap” Funktion beim Mac kann den RCommander empfindlich ausbremsen. Schalten Sie diese Funktion aus z.B. im RCommander über Tools - Manage Mac OS X app nap for R.app.

Sei aktuell

  • Verwenden Sie möglichst die neueste Version von R, RStudio und Ihres Betriebssystems. Ältere Versionen führen u.U. zu Problemen; je älter, desto Problem…

  • Updaten Sie Ihre Packages regelmäßig z.B. mit update.packages() oder dem Button “Update” bei RStudio (Reiter Packages).

Sei milde

R zu lernen kann hart sein. Ich weiß, wovon ich spreche. Wahrscheinlich eine spirituelle Prüfung in Geduld und Hartnäckigkeit… Tolle Gelegenheit, sich in diesen Tugenden zu trainieren :-)

Thanks to my student Marie Halbich who took the pains to collect the data!

At times, your data set will be in “wide” format, i.e, many columns in comparison to rows. For some analyses however, it is more suitable to have the data in “long” format. That is, many rows in comparison to columns.

Let’s have a look at this data set, for example.

d <- read.csv("https://sebastiansauer.github.io/data/facial_beauty_raw.csv")

This is the data from a study tapping into the effect of computerized “beautification” of some faces on subjective “like”.

In short, two strawpersons were present either in “natural form” or “beautified” with some computer help. Note that the design is a between-group design with the factor “strawpersons” (A and B) and the factor “presentation” (natural vs. beautified).

library(dplyr)
glimpse(d)
#> Observations: 1,440
#> Variables: 23
#> $ X                                          <int> 1, 2, 3, 4, 5, 6, 7...
#> $ Item                                       <int> 1, 2, 3, 4, 5, 6, 7...
#> $ Alter                                      <int> 26, 31, 35, 34, 31,...
#> $ girl_C_unbearbeitet_Likes                  <int> NA, 1, 1, 1, 1, 1, ...
#> $ girl_C_unbearbeitet_Dislikes               <int> 1, NA, NA, NA, NA, ...
#> $ girl_C_unbearbeitet_Superlikes             <int> NA, NA, NA, NA, NA,...
#> $ girl_C_unbearbeitet_Kontaktaufnahmen       <int> NA, 1, 1, NA, 1, 1,...
#> $ girl_C_unbearbeitet_keine_Kontaktaufnahmen <int> 1, NA, NA, 1, NA, N...
#> $ girl_C_bearbeitet_Like                     <int> NA, NA, NA, NA, NA,...
#> $ girl_C_bearbeitet_Dislike                  <int> NA, NA, NA, NA, NA,...
#> $ girl_C_bearbeitet_Superlike                <int> NA, NA, NA, NA, NA,...
#> $ girl_C_bearbeitet_Kontaktaufnahme          <int> NA, NA, NA, NA, NA,...
#> $ girl_C_bearbeitet_keine_Kontaktaufnahme    <int> NA, NA, NA, NA, NA,...
#> $ girl_A_unbearbeitet_Like                   <int> NA, NA, NA, NA, NA,...
#> $ girl_A_unbearbeitet_Dislike                <int> NA, NA, NA, NA, NA,...
#> $ girl_A_unbearbeitet_Superlike              <int> NA, NA, NA, NA, NA,...
#> $ girl_A_unbearbeitet_Kontaktaufnahme        <int> NA, NA, NA, NA, NA,...
#> $ girl_A_unbearbeitet_keine_Kontaktaufnahme  <int> NA, NA, NA, NA, NA,...
#> $ girl_A_bearbeitet_Like                     <int> NA, NA, NA, NA, NA,...
#> $ girl_A_bearbeitet_Dislike                  <int> NA, NA, NA, NA, NA,...
#> $ girl_A_bearbeitet_Superlike                <int> NA, NA, NA, NA, NA,...
#> $ girl_A_bearbeitet_Kontaktaufnahme          <int> NA, NA, NA, NA, NA,...
#> $ girl_A_bearbeitet_keine_Kontaktaufnahme    <int> NA, NA, NA, NA, NA,...

On closer inspection, we recognize a great number of missing values. In fact, the data frame is structured like this:

library(knitr)
include_graphics("facial_beauty.png")

plot of chunk unnamed-chunk-3

Where each blue rectangle represents the “core” data set for one of the four conditions (mentioned above). All the grey area represents vast deserts of NAs.

That said, the data set would be much nicer if the four “core data sets” would be aligned beneath each other like this:

library(knitr)
include_graphics("beauty_aligned.png")

plot of chunk unnamed-chunk-4

Of course, this would demand that the columns be the same in each blue square. And we will need a column indicating the value of the factor “strawperson” and a columnd for “presentation”. We don’t need column headers for times, only once, as shown on the diagram.

Cutting out the sub data frames

So we will “cut out” each sub data frame (blue rectangle) and stick them together one beneath another.

Let’s start for “girl C” and “unbearbeitet” (meaning “raw” or “natural”).

d_C_raw <- 
  d %>% 
  dplyr::select(dplyr::contains("girl_C_unbearbeitet"), Alter)

The helper function contains is handy to choose the right columns.

Now we get rid of the empty lines. We exclude all rows where only NAs appear (not looking at Alter, were no NAs will pop up).

d_C_raw <- 
  d_C_raw %>% 
    filter(!is.na(girl_C_unbearbeitet_Likes) |
          !is.na(girl_C_unbearbeitet_Dislikes) |
          !is.na(girl_C_unbearbeitet_Superlikes)| 
            !is.na(girl_C_unbearbeitet_Kontaktaufnahmen) |
            !is.na(girl_C_unbearbeitet_keine_Kontaktaufnahmen))

Note that the | operator means OR (logical OR). So we way that we want any raw where there is no NA at girl_C_unbearbeitet_Likes or no NA at girl_C_unbearbeitet_Dislikes or no NA at girl_C_unbearbeitet_Superlikes.

Now we have our first “data cubicle”. Let’s repeat that for the other 3 data cubicles.

Data frame for girl C “beautified”:

d_C_beauty <-
  d %>% 
  dplyr::select(dplyr::contains("girl_C_bearbeitet"), Alter)

d_C_beauty <- 
 d_C_beauty %>% 
    filter(!is.na(girl_C_bearbeitet_Like) |
          !is.na(girl_C_bearbeitet_Dislike) |
          !is.na(girl_C_bearbeitet_Superlike) | 
            !is.na(girl_C_bearbeitet_Kontaktaufnahme) |
            !is.na(girl_C_bearbeitet_keine_Kontaktaufnahme))

The variable names are probably somewhat verbose, but let’s not discuss that here. Also note that the names are not fully consistent; at times they come in plural and at times in singular form.

Now girl A, natural:

d_A_raw <-
  d %>% 
  dplyr::select(dplyr::contains("girl_A_unbearbeitet"), Alter)

d_A_raw <- 
 d_A_raw %>% 
    filter(!is.na(girl_A_unbearbeitet_Like) |
          !is.na(girl_A_unbearbeitet_Dislike) |
          !is.na(girl_A_unbearbeitet_Superlike) | 
            !is.na(girl_A_unbearbeitet_Kontaktaufnahme) |
            !is.na(girl_A_unbearbeitet_keine_Kontaktaufnahme))

Once more for girl A, processed:

d_A_beauty <-
  d %>% 
  dplyr::select(dplyr::contains("girl_A_bearbeitet"), Alter)

d_A_beauty <- 
 d_A_beauty %>% 
    filter(!is.na(girl_A_bearbeitet_Like) |
          !is.na(girl_A_bearbeitet_Dislike) |
          !is.na(girl_A_bearbeitet_Superlike) | 
            !is.na(girl_A_bearbeitet_Kontaktaufnahme) |
            !is.na(girl_A_bearbeitet_keine_Kontaktaufnahme))

Adjust names of sub data frames

Before we can bind the sub data frames together, we have to make sure the names of the columns are identical. So let’s do that now. In addition, we need to save the information about the study factors (natural vs. processed; girl A vs. girl C).

d_A_beauty$girl <- "A"
d_A_beauty$presentation <- "beauty"
d_A_raw$girl <- "A"
d_A_raw$presentation <- "natural"


d_C_beauty$girl <- "C"
d_C_beauty$presentation <- "beauty"
d_C_raw$girl <- "C"
d_C_raw$presentation <- "natural"

OK, now the names of the variables:

names_cols <- c("like_yes", "like_no", "superlike_yes", "contact_yes", "contact_no", "age", "girl", "presentation")

names(d_A_raw) <- names_cols
names(d_A_beauty) <- names_cols
names(d_C_raw) <- names_cols
names(d_C_beauty) <- names_cols

Bind sub data frames

Now we are ready to bind the four sub data frames into one, by placing the four pieces underneath each other:

d_long <- bind_rows(d_A_raw, d_A_beauty, d_C_raw, d_C_beauty)

Wow! Now we have reformated our data frame from “wide” format to “long” format. Long ride.

Unite columns

Note that some columns are dependent: a 1 in like corresponds to NA in dislike and a NA in superlike. We can merge these columns into one. Same story for contact.

First, recode NAs to 0 as this it what they mean here.


d_long <-
  d_long %>% 
  mutate_at(vars(like_yes, like_no, superlike_yes, contact_yes, contact_no), 
            recode, .missing = 0, `1` = 1)

A quick check:

count(d_long, like_yes)
#> # A tibble: 2 × 2
#>   like_yes     n
#>      <dbl> <int>
#> 1        0   404
#> 2        1  1036
count(d_long,like_no)
#> # A tibble: 2 × 2
#>   like_no     n
#>     <dbl> <int>
#> 1       0  1036
#> 2       1   404

Seems OK.

Next, merge the three “like columns” into one using dplyr::case_when:


d_long <-
  d_long %>% 
  mutate(like = case_when(
    d_long$like_yes == 1 ~ "like",
    d_long$like_no == 1 ~ "dislike",
    d_long$superlike_yes == 1 ~ "superlike"))

Note that case_when will fail if there are NAs in your vector (without an helpful error).

Now we can drop the initial three “like columns”.

d_long <- 
  d_long %>% 
  select(-c(like_yes, like_no, superlike_yes))

And now same story for the two “contact columns”.

Wait, if 1 in contact_yes means 0 in contact_no, then we are already done. Let’s look at the counts of the two variables:

count(d_long, contact_yes)
#> # A tibble: 2 × 2
#>   contact_yes     n
#>         <dbl> <int>
#> 1           0   853
#> 2           1   587
count(d_long, contact_no)
#> # A tibble: 2 × 2
#>   contact_no     n
#>        <dbl> <int>
#> 1          0   585
#> 2          1   855

Hm, that’s no exact match, but nearly. Probably some typing errors or something similar. So we can just rename one variable, drop the other, and that’s it.

d_long <-
  d_long %>% 
  rename(contact = contact_yes) %>% 
  select(-contact_no)

Debrief

That’s about it. We have now a data frame in “long” format. This version is more adequat for many analyses. Have fun!

Session info

sessionInfo()
#> R version 3.3.2 (2016-10-31)
#> Platform: x86_64-apple-darwin13.4.0 (64-bit)
#> Running under: macOS Sierra 10.12.2
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] knitr_1.15.1 dplyr_0.5.0 
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_0.12.8         codetools_0.2-15    digest_0.6.11      
#>  [4] rprojroot_1.1       assertthat_0.1      R6_2.2.0           
#>  [7] DBI_0.5-1           backports_1.0.4     magrittr_1.5       
#> [10] evaluate_0.10       highr_0.6           stringi_1.1.2      
#> [13] lazyeval_0.2.0.9000 rmarkdown_1.3       tools_3.3.2        
#> [16] stringr_1.1.0       rsconnect_0.7       yaml_2.1.14        
#> [19] htmltools_0.3.5     tibble_1.2


This YACSDA (Yet-another-case-study-on-data-analysis) in composed in German language. Some typical data analytical steps are introduced.


Wovon ist die Häufigkeit von Affären (Seitensprüngen) in Ehen abhängig? Diese Frage soll anhand des Datensates Affair untersucht werden.

Dieser Post stellt beispielhaft eine grundlegende Methoden der praktischen Datenanalyse im Rahmen einer kleinen Fallstudie (YACSDA) vor.

Quelle der Daten: http://statsmodels.sourceforge.net/0.5.0/datasets/generated/fair.html

Der Datensatz findet sich (in ähnlicher Form) auch im R-Paket COUNT (https://cran.r-project.org/web/packages/COUNT/index.html).

Laden wir als erstes den Datensatz in R. Wählen Sie zuerst das Verzeichnis als Arbeitsverzeichnis, in dem die Daten liegen. Dann laden Sie z.B. mit dem R-Commander (s. Skript) oder “per Hand” z.B. bei mir so:

Affair <- read.csv("~/Documents/Literatur/Methoden_Literatur/Datensaetze/Affair.csv")

Schauen wir mal, ob es funktioniert hat (“Datenmatrix betrachten”):

head(Affair)
#>   affairs gender age yearsmarried children religiousness education
#> 1       0   male  37        10.00       no             3        18
#> 2       0 female  27         4.00       no             4        14
#> 3       0 female  32        15.00      yes             1        12
#> 4       0   male  57        15.00      yes             5        18
#> 5       0   male  22         0.75       no             2        17
#> 6       0 female  32         1.50       no             2        17
#>   occupation rating
#> 1          7      4
#> 2          6      4
#> 3          1      4
#> 4          6      5
#> 5          6      3
#> 6          5      5

Ok scheint zu passen. Was jetzt?

Geben Sie zentrale deskriptive Statistiken an für Affärenhäufigkeit und Ehezufriedenheit!

# nicht robust:
mean(Affair$affairs, na.rm = T)
#> [1] 1.455907
sd(Affair$affairs, na.rm = T)
#> [1] 3.298758
# robust:
median(Affair$affair, na.rm = T)
#> [1] 0
IQR(Affair$affair, na.rm = T)
#> [1] 0

Es scheint, die meisten Leute haben keine Affären:

table(Affair$affairs)
#> 
#>   0   1   2   3   7  12 
#> 451  34  17  19  42  38

Man kann sich viele Statistiken mit dem Befehl describe aus psych ausgeben lassen, das ist etwas praktischer:

library(psych)
                 
describe(Affair$affairs)
#>    vars   n mean  sd median trimmed mad min max range skew kurtosis   se
#> X1    1 601 1.46 3.3      0    0.55   0   0  12    12 2.34     4.19 0.13
describe(Affair$rating)
#>    vars   n mean  sd median trimmed  mad min max range  skew kurtosis   se
#> X1    1 601 3.93 1.1      4    4.07 1.48   1   5     4 -0.83    -0.22 0.04

Dazu muss das Paket psych natürlich vorher installiert sein. Beachten Sie, dass man ein Paket nur einmal installieren muss, aber jedes Mal, wen Sie R starten, auch starten muss (mit library).

install.packages("psych")

Visualisieren Sie zentrale Variablen!

Sicherlich sind Diagramme auch hilfreich. Dies geht wiederum mit dem R-Commander oder z.B. mit folgenden Befehlen:


library(ggplot2)
qplot(x = affairs, data = Affair)
qplot(x = rating, data = Affair)

plot of chunk unnamed-chunk-7plot of chunk unnamed-chunk-7

Die meisten Menschen (dieser Stichprobe) scheinen mit Ihrer Beziehung sehr zufrieden zu sein.

Wer ist zufriedener mit der Partnerschaft: Personen mit Kindern oder ohne?

Nehmen wir dazu mal ein paar dplyr-Befehle:

library(dplyr)
Affair %>% 
  group_by(children) %>% 
  summarise(rating_children = mean(rating, na.rm = T))
#> # A tibble: 2 × 2
#>   children rating_children
#>     <fctr>           <dbl>
#> 1       no        4.274854
#> 2      yes        3.795349

Ah! Kinder sind also ein Risikofaktor für eine Partnerschaft! Gut, dass wir das geklärt haben.

Wie viele fehlende Werte gibt es? Was machen wir am besten damit?

Diesen Befehl könnten wir für jede Spalte auführen:

sum(is.na(Affair$affairs))
#> [1] 0

Oder lieber alle auf einmal:

Affair %>% 
  summarise_each(funs(sum(is.na(.))))
#>   affairs gender age yearsmarried children religiousness education
#> 1       0      0   0            0        0             0         0
#>   occupation rating
#> 1          0      0

Übrigens gibt es ein gutes Cheat Sheet für dplyr.

Ah, gut, keine fehlenden Werte. Das macht uns das Leben leichter.

Wer ist glücklicher: Männer oder Frauen?

Affair %>% 
  group_by(gender) %>% 
  summarise(rating_gender = mean(rating))
#> # A tibble: 2 × 2
#>   gender rating_gender
#>   <fctr>         <dbl>
#> 1 female      3.939683
#> 2   male      3.923077

Praktisch kein Unterschied. Heißt das auch, es gibt keinen Unterschied in der Häufigkeit der Affären?

Affair %>% 
  group_by(gender) %>% 
  summarise(affairs_gender = mean(affairs))
#> # A tibble: 2 × 2
#>   gender affairs_gender
#>   <fctr>          <dbl>
#> 1 female       1.419048
#> 2   male       1.496503

Scheint auch kein Unterschied zu sein…

Und zum Abschluss noch mal etwas genauer: Teilen wir mal nach Geschlecht und nach Kinderstatus auf, also in 4 Gruppen. Theoretisch dürfte es hier auch keine Unterschiede/Zusammenhänge geben. Zumindest fällt mir kein sinnvoller Grund ein; zumal die vorherige eindimensionale Analyse keine Unterschiede zu Tage gefördert hat.

Affair %>% 
  group_by(gender, children) %>% 
  summarise(affairs_mean = mean(affairs),
            rating_mean = mean(rating))
#> Source: local data frame [4 x 4]
#> Groups: gender [?]
#> 
#>   gender children affairs_mean rating_mean
#>   <fctr>   <fctr>        <dbl>       <dbl>
#> 1 female       no    0.8383838    4.404040
#> 2 female      yes    1.6851852    3.726852
#> 3   male       no    1.0138889    4.097222
#> 4   male      yes    1.6588785    3.864486

Affair %>% 
  group_by(children, gender) %>% 
  summarise(affairs_mean = mean(affairs),
            rating_mean = mean(rating))
#> Source: local data frame [4 x 4]
#> Groups: children [?]
#> 
#>   children gender affairs_mean rating_mean
#>     <fctr> <fctr>        <dbl>       <dbl>
#> 1       no female    0.8383838    4.404040
#> 2       no   male    1.0138889    4.097222
#> 3      yes female    1.6851852    3.726852
#> 4      yes   male    1.6588785    3.864486

Berichten Sie eine relevante Effektstärke!

Hm, auch keine gewaltigen Unterschiede. Höchstens für die Zufriedenheit mit der Partnerschaft bei kinderlosen Personen scheinen sich Männer und Frauen etwas zu unterscheiden. Hier stellt sich die Frage nach der Größe des Effekts, z.B. anhand Cohen’s d. Dafür müssen wir noch die SD pro Gruppe wissen:

Affair %>% 
  group_by(children, gender) %>% 
  summarise(rating_mean = mean(rating),
            rating_sd = sd(rating))
#> Source: local data frame [4 x 4]
#> Groups: children [?]
#> 
#>   children gender rating_mean rating_sd
#>     <fctr> <fctr>       <dbl>     <dbl>
#> 1       no female    4.404040 0.9138302
#> 2       no   male    4.097222 1.0636070
#> 3      yes female    3.726852 1.1829884
#> 4      yes   male    3.864486 1.0460525
d <- (4.4 - 4.1)/(1)

Die Effektstärke beträgt etwa 0.3.

Berechnen und visualisieren Sie zentrale Korrelationen!

Affair %>% 
  select_if(is.numeric) %>% 
  cor -> cor_tab

cor_tab
#>                    affairs        age yearsmarried religiousness
#> affairs        1.000000000  0.0952372   0.18684169   -0.14450135
#> age            0.095237204  1.0000000   0.77754585    0.19377693
#> yearsmarried   0.186841686  0.7775458   1.00000000    0.21826067
#> religiousness -0.144501345  0.1937769   0.21826067    1.00000000
#> education     -0.002437441  0.1345960   0.04000272   -0.04257108
#> occupation     0.049611758  0.1664125   0.04459201   -0.03972232
#> rating        -0.279512403 -0.1989999  -0.24311883    0.02429578
#>                  education  occupation      rating
#> affairs       -0.002437441  0.04961176 -0.27951240
#> age            0.134596015  0.16641254 -0.19899990
#> yearsmarried   0.040002716  0.04459201 -0.24311883
#> religiousness -0.042571079 -0.03972232  0.02429578
#> education      1.000000000  0.53360524  0.10930347
#> occupation     0.533605242  1.00000000  0.01742227
#> rating         0.109303473  0.01742227  1.00000000

library(corrplot)
corrplot(cor_tab)

plot of chunk unnamed-chunk-16

Wie groß ist der Einfluss (das Einflussgewicht) der Ehejahre bzw. Ehezufriedenheit auf die Anzahl der Affären?

Dazu sagen wir R: “Hey R, rechne mal ein lineares Modell”, also eine normale (lineare) Regression. Dazu können wir entweder das entsprechende Menü im R-Commander auswählen, oder folgende R-Befehle ausführen:

lm1 <- lm(affairs ~ yearsmarried, data = Affair)
summary(lm1)  # Ergebnisse der Regression zeigen
#> 
#> Call:
#> lm(formula = affairs ~ yearsmarried, data = Affair)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2.2106 -1.6575 -0.9937 -0.5974 11.3658 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   0.55122    0.23511   2.345   0.0194 *  
#> yearsmarried  0.11063    0.02377   4.655    4e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.243 on 599 degrees of freedom
#> Multiple R-squared:  0.03491,	Adjusted R-squared:  0.0333 
#> F-statistic: 21.67 on 1 and 599 DF,  p-value: 3.996e-06
lm2 <- lm(affairs ~ rating, data = Affair)
summary(lm2)
#> 
#> Call:
#> lm(formula = affairs ~ rating, data = Affair)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.9063 -1.3989 -0.5631 -0.5631 11.4369 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   4.7421     0.4790   9.900   <2e-16 ***
#> rating       -0.8358     0.1173  -7.125    3e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.17 on 599 degrees of freedom
#> Multiple R-squared:  0.07813,	Adjusted R-squared:  0.07659 
#> F-statistic: 50.76 on 1 and 599 DF,  p-value: 3.002e-12

Also: yearsmarried und rating sind beide statistisch signifikante Prädiktoren für die Häufigkeit von Affären. Das adjustierte ist allerdings in beiden Fällen nicht so groß.

Um wie viel erhöht sich die erklärte Varianz (R-Quadrat) von Affärenhäufigkeit wenn man den Prädiktor Ehezufriedenheit zum Prädiktor Ehejahre hinzufügt? (Wie) verändern sich die Einflussgewichte (b)?

lm3 <- lm(affairs ~ rating + yearsmarried, data = Affair)
lm4 <- lm(affairs ~ yearsmarried + rating, data = Affair)
summary(lm3)
#> 
#> Call:
#> lm(formula = affairs ~ rating + yearsmarried, data = Affair)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4.1474 -1.6495 -0.8365 -0.1616 11.8945 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   3.76913    0.56715   6.646 6.80e-11 ***
#> rating       -0.74395    0.12005  -6.197 1.07e-09 ***
#> yearsmarried  0.07481    0.02377   3.147  0.00173 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.147 on 598 degrees of freedom
#> Multiple R-squared:  0.09315,	Adjusted R-squared:  0.09012 
#> F-statistic: 30.71 on 2 and 598 DF,  p-value: 2.01e-13
summary(lm4)
#> 
#> Call:
#> lm(formula = affairs ~ yearsmarried + rating, data = Affair)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4.1474 -1.6495 -0.8365 -0.1616 11.8945 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   3.76913    0.56715   6.646 6.80e-11 ***
#> yearsmarried  0.07481    0.02377   3.147  0.00173 ** 
#> rating       -0.74395    0.12005  -6.197 1.07e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.147 on 598 degrees of freedom
#> Multiple R-squared:  0.09315,	Adjusted R-squared:  0.09012 
#> F-statistic: 30.71 on 2 and 598 DF,  p-value: 2.01e-13

Ok. Macht eigentlich die Reihenfolge der Prädiktoren in der Regression einen Unterschied? Der Vergleich von Modell 3 vs. Modell 4 beantwortet diese Frage.

Wir sehen, dass beim 1. Regressionsmodell das R^2 0.03 war; beim 2. Modell 0.08 und beim 3. Modell liegt R^2 bei 0.09. Die Differenz zwischen Modell 1 und 3 liegt bei (gerundet) 0.06; wenig.

Welche Prädiktoren würden Sie noch in die Regressionsanalyse aufnehmen?

Hm, diese Frage klingt nicht so, als ob der Dozent die Antwort selber wüsste… Naja, welche Variablen gibt es denn alles:

#> [1] "affairs"       "gender"        "age"           "yearsmarried" 
#> [5] "children"      "religiousness" "education"     "occupation"   
#> [9] "rating"

Z.B. wäre doch interessant, ob Ehen mit Kinder mehr oder weniger Seitensprüngen aufweisen. Und ob die “Kinderfrage” die anderen Zusammenhänge/Einflussgewichte in der Regression verändert. Probieren wir es auch. Wir können wiederum im R-Comamnder ein Regressionsmodell anfordern oder es mit der Syntax probieren:

lm5 <- lm(affairs~ rating + yearsmarried + children, data = Affair)
summary(lm5)
#> 
#> Call:
#> lm(formula = affairs ~ rating + yearsmarried + children, data = Affair)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4.3537 -1.7316 -0.8927 -0.1719 12.0162 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   3.85245    0.58808   6.551 1.24e-10 ***
#> rating       -0.74861    0.12043  -6.216 9.57e-10 ***
#> yearsmarried  0.08332    0.02853   2.921  0.00362 ** 
#> childrenyes  -0.18805    0.34817  -0.540  0.58932    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.148 on 597 degrees of freedom
#> Multiple R-squared:  0.09359,	Adjusted R-squared:  0.08904 
#> F-statistic: 20.55 on 3 and 597 DF,  p-value: 1.107e-12
r2_lm5 <- summary(lm5)$r.squared

Das Regressionsgewicht von childrenyes ist negativ. Das bedeutet, dass Ehen mit Kindern weniger Affären verbuchen (aber geringe Zufriedenheit, wie wir oben gesehen haben! Hrks!). Allerdings ist der p-Wert nich signifikant, was wir als Zeichen der Unbedeutsamkeit dieses Prädiktors verstehen können. lungert immer noch bei mickrigen 0.09 herum. Wir haben bisher kaum verstanden, wie es zu Affären kommt. Oder unsere Daten bergen diese Informationen einfach nicht.

Wir könnten auch einfach mal Prädiktoren, die wir haben, ins Feld schicken. Mal sehen, was dann passiert:

lm6 <- lm(affairs ~ ., data = Affair)
summary(lm6)
#> 
#> Call:
#> lm(formula = affairs ~ ., data = Affair)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -5.0503 -1.7226 -0.7947  0.2101 12.7036 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)    5.87201    1.13750   5.162 3.34e-07 ***
#> gendermale     0.05409    0.30049   0.180   0.8572    
#> age           -0.05098    0.02262  -2.254   0.0246 *  
#> yearsmarried   0.16947    0.04122   4.111 4.50e-05 ***
#> childrenyes   -0.14262    0.35020  -0.407   0.6840    
#> religiousness -0.47761    0.11173  -4.275 2.23e-05 ***
#> education     -0.01375    0.06414  -0.214   0.8303    
#> occupation     0.10492    0.08888   1.180   0.2383    
#> rating        -0.71188    0.12001  -5.932 5.09e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.095 on 592 degrees of freedom
#> Multiple R-squared:  0.1317,	Adjusted R-squared:   0.12 
#> F-statistic: 11.23 on 8 and 592 DF,  p-value: 7.472e-15
r2_lm6 <- round(summary(lm6)$r.squared, 2)

Der “.” im Befehl affairs ~ . oben soll sagen: nimm “alle Variablen, die noch in der Datenmatrix übrig sind”.

Insgesamt bleibt die erklärte Varian in sehr bescheidenem Rahmen: 0.13. Das zeigt uns, dass es immer noch nur schlecht verstanden ist – im Rahmen dieser Analyse – welche Faktoren die Affärenhäufigkeit erklärt.

Unterscheiden sich die Geschlechter statistisch signifikant? Wie groß ist der Unterschied? Sollte hierlieber das d-Maß oder Rohwerte als Effektmaß angegeben werden?

Hier bietet sich ein t-Test für unabhängige Gruppen an. Die Frage lässt auf eine ungerichtete Hypothese schließen ( sei .05). Mit dem entsprechenden Menüpunkt im R-Commander oder mit folgender Syntax lässt sich diese Analyse angehen:

t1 <- t.test(affairs ~ gender, data = Affair)
t1
#> 
#> 	Welch Two Sample t-test
#> 
#> data:  affairs by gender
#> t = -0.28733, df = 594.01, p-value = 0.774
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -0.6068861  0.4519744
#> sample estimates:
#> mean in group female   mean in group male 
#>             1.419048             1.496503

Der p-Wert ist mit 0.7739606 > . Daher wird die beibehalten. Auf Basis der Stichprobendaten entscheiden wir uns für die . Entsprechend umschließt das 95%-KI die Null.

Da die Differenz nicht signifikant ist, kann argumentiert werden, dass wir d auf 0 schätzen müssen. Man kann sich den d-Wert auch z.B. von {MBESS} schätzen lassen.

Dafür brauchen wir die Anzahl an Männer und Frauen: 315, 286.

library(MBESS)
ci.smd(ncp = t1$statistic,
    n.1 = 315,
    n.2 = 286)
#> $Lower.Conf.Limit.smd
#> [1] -0.1835475
#> 
#> $smd
#>           t 
#> -0.02346813 
#> 
#> $Upper.Conf.Limit.smd
#> [1] 0.1366308

Das Konfidenzintervall ist zwar relativ klein (die Schätzung also aufgrund der recht großen Stichprobe relativ präzise), aber der Schätzwert für d smd liegt sehr nahe bei Null. Das stärkt unsere Entscheidung, von einer Gleichheit der Populationen (Männer vs. Frauen) auszugehen.

Rechnen Sie die Regressionsanalyse getrennt für kinderlose Ehe und Ehen mit Kindern!

Hier geht es im ersten Schritt darum, die entsprechenden Teil-Mengen der Datenmatrix zu erstellen. Das kann man natürlich mit Excel o.ä. tun. Alternativ könnte man es in R z.B. so machen:

Affair2 <- Affair[Affair$children == "yes", ]
lm7 <- lm(affairs~ rating, data = Affair2)
summary(lm7)
#> 
#> Call:
#> lm(formula = affairs ~ rating, data = Affair2)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4.1903 -1.4877 -0.5869 -0.4877 11.4131 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   5.0912     0.5701   8.930  < 2e-16 ***
#> rating       -0.9009     0.1441  -6.252 9.84e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.336 on 428 degrees of freedom
#> Multiple R-squared:  0.08367,	Adjusted R-squared:  0.08153 
#> F-statistic: 39.08 on 1 and 428 DF,  p-value: 9.845e-10

Affair3 <- Affair[Affair$children == "no", ]
lm8 <- lm(affairs~ rating, data = Affair3)
summary(lm8)
#> 
#> Call:
#> lm(formula = affairs ~ rating, data = Affair3)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2.5465 -1.0494 -0.5504 -0.5504 11.4496 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)   
#> (Intercept)   3.0455     0.9140   3.332  0.00106 **
#> rating       -0.4990     0.2083  -2.395  0.01771 * 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2.685 on 169 degrees of freedom
#> Multiple R-squared:  0.03283,	Adjusted R-squared:  0.02711 
#> F-statistic: 5.737 on 1 and 169 DF,  p-value: 0.01771

Übrigens, einfacher geht das “Subsetten” so:

library(dplyr)
Affair4 <- filter(Affair, children == "yes")
head(Affair4)
#>   affairs gender age yearsmarried children religiousness education
#> 1       0 female  32           15      yes             1        12
#> 2       0   male  57           15      yes             5        18
#> 3       0   male  57           15      yes             2        14
#> 4       0 female  32           15      yes             4        16
#> 5       0   male  37           15      yes             2        20
#> 6       0   male  27            4      yes             4        18
#>   occupation rating
#> 1          1      4
#> 2          6      5
#> 3          4      4
#> 4          1      2
#> 5          7      2
#> 6          6      4

Rechnen Sie die Regression nur für “Halodries”; d.h. für Menschen mit Seitensprüngen. Dafür müssen Sie alle Menschen ohne Affären aus den Datensatz entfernen.

Also, rechnen wir nochmal die Standardregression (lm1). Probieren wir den Befehl filter dazu nochmal aus:

Affair5 <- filter(Affair, affairs != 0)
lm9 <- lm(affairs ~ rating, data = Affair5)
summary(lm9)
#> 
#> Call:
#> lm(formula = affairs ~ rating, data = Affair5)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.0605 -3.5157 -0.0605  3.6895  7.4843 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   8.7570     1.0225   8.564  1.3e-14 ***
#> rating       -0.8483     0.2800  -3.030  0.00289 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 4.144 on 148 degrees of freedom
#> Multiple R-squared:  0.05841,	Adjusted R-squared:  0.05205 
#> F-statistic: 9.181 on 1 and 148 DF,  p-value: 0.002887

Berechnen Sie für eine logistische Regression mit “Affäre ja vs. nein” als Kriterium, wie stark der Einfluss von Geschlecht, Kinderstatus, Ehezufriedenheit und Ehedauer ist!


Affair %>% 
  mutate(affairs_dichotom = if_else(affairs == 0, 0, 1)) %>% 
  glm(affairs_dichotom ~gender + children + rating + yearsmarried, data = .) -> lm10

summary(lm10)
#> 
#> Call:
#> glm(formula = affairs_dichotom ~ gender + children + rating + 
#>     yearsmarried, data = .)
#> 
#> Deviance Residuals: 
#>      Min        1Q    Median        3Q       Max  
#> -0.57662  -0.26766  -0.17186  -0.06459   0.93295  
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   0.515778   0.079344   6.501 1.69e-10 ***
#> gendermale    0.037945   0.034218   1.109    0.268    
#> childrenyes   0.054032   0.046307   1.167    0.244    
#> rating       -0.090337   0.015985  -5.651 2.47e-08 ***
#> yearsmarried  0.003947   0.003787   1.042    0.298    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 0.1746456)
#> 
#>     Null deviance: 112.56  on 600  degrees of freedom
#> Residual deviance: 104.09  on 596  degrees of freedom
#> AIC: 663.8
#> 
#> Number of Fisher Scoring iterations: 2

Wenn if_else unbekannt ist, lohnt sich ein Blick in die Hilfe mit ?if_else (dplyr muss vorher geladen sein).

Aha, signifikant ist die Ehezufriedenheit: Je größer rating desto geringer die Wahrscheinlickeit für affairs_dichotom. Macht Sinn! Je zufriedener, desto geringer Seitensprung-Gefahr…

Übrigens, die Funktion lm und glm spucken leider keine brave Tabelle in Normalform aus. Aber man leicht eine schöne Tabelle (data.frame) bekommen mit dem Befehl tidy aus broom:

library(broom)
tidy(lm10) 
#>           term     estimate   std.error statistic      p.value
#> 1  (Intercept)  0.515777667 0.079343849  6.500537 1.693109e-10
#> 2   gendermale  0.037944509 0.034217912  1.108908 2.679172e-01
#> 3  childrenyes  0.054032156 0.046307172  1.166820 2.437495e-01
#> 4       rating -0.090336904 0.015984908 -5.651387 2.468646e-08
#> 5 yearsmarried  0.003946811 0.003786673  1.042290 2.976998e-01

Und Tabellen (d.h. brave Dataframes) kann man sich schön ausgeben lassen z.B. mit dem Befehl knitr::kable:

library(knitr)
tidy(lm10) %>% kable
term estimate std.error statistic p.value
(Intercept) 0.5157777 0.0793438 6.500538 0.0000000
gendermale 0.0379445 0.0342179 1.108908 0.2679172
childrenyes 0.0540322 0.0463072 1.166821 0.2437495
rating -0.0903369 0.0159849 -5.651387 0.0000000
yearsmarried 0.0039468 0.0037867 1.042290 0.2976998

Visualisieren wir mal was!

Ok, wie wäre es damit:

Affair %>% 
   select(affairs, gender, children, rating) %>%
  ggplot(aes(x = affairs, y = rating)) + geom_jitter(aes(color = gender, shape = children)) 

plot of chunk unnamed-chunk-31

Affair %>% 
   mutate(rating_dichotom = ntile(rating, 2)) %>% 
   ggplot(aes(x = yearsmarried, y = affairs)) + geom_jitter(aes(color = gender)) +
  geom_smooth()

plot of chunk unnamed-chunk-32

Puh. Geschafft!

Versionshinweise und SessionInfo

  • Datum erstellt: 2017-01-06
  • R Version: 3.3.2
  • dplyr Version: 0.5.0
sessionInfo()
#> R version 3.3.2 (2016-10-31)
#> Platform: x86_64-apple-darwin13.4.0 (64-bit)
#> Running under: macOS Sierra 10.12.2
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] broom_0.4.1   MBESS_4.1.0   corrplot_0.77 ggplot2_2.2.1 psych_1.6.9  
#> [6] knitr_1.15.1  dplyr_0.5.0  
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_0.12.8         magrittr_1.5        munsell_0.4.3      
#>  [4] mnormt_1.5-5        lattice_0.20-34     colorspace_1.3-2   
#>  [7] R6_2.2.0            plyr_1.8.4          stringr_1.1.0      
#> [10] highr_0.6           tools_3.3.2         parallel_3.3.2     
#> [13] grid_3.3.2          nlme_3.1-128        gtable_0.2.0       
#> [16] DBI_0.5-1           htmltools_0.3.5     yaml_2.1.14        
#> [19] lazyeval_0.2.0.9000 assertthat_0.1      rprojroot_1.1      
#> [22] digest_0.6.11       tibble_1.2          tidyr_0.6.0        
#> [25] reshape2_1.4.2      codetools_0.2-15    rsconnect_0.7      
#> [28] evaluate_0.10       rmarkdown_1.3       labeling_0.3       
#> [31] stringi_1.1.2       scales_0.4.1        backports_1.0.4    
#> [34] foreign_0.8-67

The variance of some data can be defined in rough terms as the mean of the squared deviations from the mean.

Let’s repeat that because it is important:

Variance: Mean of squared deviations from the mean.

An example helps to illustrate. Assume some class of students are forced to write an exam in a statistics class (OMG). Let’s say the grades range fom 1 to 6, 1 being the best and 6 the worst. We can compute the mean of this class (eg., 2.3); once we know the mean, we can subtract the mean from each individual grade. If Anna scored a 3 (OK but not exciting), “her” residual, delta or deviation would 3 - 2.3 = 0.7, and so on. The following picture illustrates this example.

library(tidyverse)

# df with students' grades (12 students)
grades_df <- data_frame(grades = c(3,2,4,1,1,2,5,1,1,3,2,2))


# compute deltas (deviations from mean)
grades_df %>% 
  mutate(grades_Delta = grades - mean(grades),
         Student = factor(row_number())) -> grades_df 

# now plot 
grades_df%>% 
  ggplot(aes(x = Student, y = grades)) + geom_point() +
  geom_hline(yintercept = mean(grades_df$grades), color = "firebrick") +
  geom_linerange(aes(x = Student, ymin = grades, ymax = grades - grades_Delta), color = "blue") 

plot of chunk unnamed-chunk-1

The vertical blue lines look a bit like little sticks…

In this example, the vertical blue sticks indicate the individual delta. Now try to envision the “average” length of the stick. That’s something similar to the variance but not exactly the same. We have to to square each stick, and then only take the avarage.

In R, a fancy (piped) way to compute the variance is this:

grades_df$grades %>% 
  `-`(mean(grades_df$grades)) %>% 
  `^`(2) %>% 
  mean
## [1] 1.520833

Literally, in plain English, this code means:

Take the data frame “grades_df”, and from there the column “grades”. Then
subtract the mean of this column. Then
take it to the 2nd power (square it, that is). Then
compute the mean of the numbers.

More formally, the variance of some empirical data is defined as

.

And for infinite variables (whatever that is) it is defined as:

.

Note that means “the mean on the long run” (infinitely long, to be sure).

Which is exactly the same as the R code above; slightly more complicated, I feel, because the steps are nested (because of the brackets). In the R code above, we “sequentialized” the steps, which renders the thing more tangible.

Note that the R function var does not divide by n but by n-1. However, for larger samples, the error is negligible. The reason is that R computes the so-called “sample variance”, that is the estimate of the population variance based on the sample data. Thus, depending on whether you are interested in a guess of the population variance or just the variance of your data at hand, the one or the other is (slightly) more appopriate. If you were to know to population variance, then again you’d divide by n(not by n-1). Here, we just divide by n and by happy.

In general, my opinion is not to worry too much about tiny details (for the purpose given), but rather to try to grasp the big pictures. So we won’t elaborate here on that further.

Note that the dice throws are independent with each others, they are not correlated.

Simulate dice throwing

With this understanding of the variance in mind, let’s continue. Let’s make R throw some dice to illustrate the additivity property of the variance.

Hey R: throw a normal die reps = 1e05 times and plot the results:

reps <- 1e05
dice <- sample(x = c(1:6), size = reps, replace = TRUE)
dice <- tibble(dice)

ggplot(dice) +
  aes(x = dice) +
  geom_histogram()

plot of chunk unnamed-chunk-3

OK, fair enough; each side appeared more or less equally often.

Now, what’ the mean and the variance?

var(dice$dice)
## [1] 2.915413
mean(dice$dice)
## [1] 3.494

OK, plausible.

Now, we are concerned with the additivity of the variance. So we are supposed to some stuff up! Suppose we repeat our dice throwing experiment, but now with 2 dices (instead of 1). After each throw, we sum up the score. After some hard thinking we feel reassured that this should yield a number between 2 and 12. Let’s have a look at the plot, as we want to get an intuitive understanding.

dice1 <- sample(x = c(1:6), size = reps, replace = TRUE)

dice2 <- sample(x = c(1:6), size = reps, replace = TRUE)

dice_sum <- dice1 + dice2
# dice_sum <- sum(dice1, dice2) this will not work!

dice_sum <- tibble(dice_sum)

ggplot(dice_sum) +
  aes(x = dice_sum) +
  geom_histogram() +
  scale_x_continuous(breaks = 1:12)  

plot of chunk unnamed-chunk-5

But what about mean and variance? Let’s see.

mean(dice_sum$dice_sum)
## [1] 7.00482
var(dice_sum$dice_sum)
## [1] 5.808515

It appears that the means are adding up, which makes sense, if you think about it. Same for the variance, it adds up!

Sum up k variables

As a last step, let’s add up, say, variables (die throws) to further strengthen our argument.

First, we make up some data with k columns of dice throws.

d <- list()
k <- 10 
for (i in 1:k){
  d[i] <- list(sample(x = c(1:6), size = reps, replace = TRUE))
}

d <- rbind.data.frame(d)
names(d) <- paste("V",1:k, sep = "")

Now sum up the score of each of the throws with k dice.

d$sum_throw <- rowSums(d[, 1:10])

We expect the variance to be about k times var(one_die_throw).

var(d$sum_throw)
## [1] 29.29128

Which is what we find. Similarly for the mean:

mean(d$sum_throw)
## [1] 34.99714

OK.

What about sd?

Does this property hold with the sd as well?

So, what’s the sd in each of the die throws?

sds <- apply(d[, 1:10], 2, sd)
sds
##       V1       V2       V3       V4       V5       V6       V7       V8 
## 1.707311 1.704684 1.707218 1.705219 1.706854 1.709697 1.709320 1.703763 
##       V9      V10 
## 1.704847 1.710907

And their sum is then 17.0698202.

But what is the sd of the sum score of the k=10 throws?

sd(d$sum_throw)
## [1] 5.412142

These two numbers are obviously not equal. So the additivity property does not hold for the sd.

More formally

Thus, we can state with some confidence:

,

where Var means variance and A, B, C, … are some variables such as dice throws.

Subtracting two variables

OK, we have seen that the variance is additive in the sense that if we sum up som arbitrary number of variables, the variance of the sum score equals the sum of the individual variances (approximately, not exactly).

But what about subtracting two variables. Does the variance get smaller in the same way, too? Let’s try and see.

dice1 <- sample(x = c(1:6), size = reps, replace = TRUE)

dice2 <- sample(x = c(1:6), size = reps, replace = TRUE)

dice_diff <- dice1 - dice2

dice_diff <- tibble(dice_diff)

ggplot(dice_diff) +
  aes(x = dice_diff) +
  geom_histogram() +
  scale_x_continuous(breaks = -5:5)  

plot of chunk unnamed-chunk-13

Obviously, if we subtract one die score from the other the result must come down between -5 and +5. But what about mean and variance? Let’s see.

mean(dice_diff$dice_diff)
## [1] -0.00232
var(dice_diff$dice_diff)
## [1] 5.813533

Interesting. The mean does what it should, it balances out again. But what about the variance? It again adds up! We see that the range is 11. Compare the “width” (range) of the individual dice throws:

ggplot(tibble(dice1)) +
  aes(x = dice1) +
  geom_histogram() +
  scale_x_continuous(breaks = 1:6)  

plot of chunk unnamed-chunk-15

The range is smaller, ie. 6. So, again, it appears that despite us subtracting variables, the variance keeps adding up!.

Formal stuff

Note that in general, the variance of some (infinite) distribution is defined as follows:

Here means “long term average” or “expected value”. X is just some variable such as a fair die throw.

This formula can be rearranged to

Mnemonic is “mean of square minus square of mean”; see here for details.

Similarly, for the sum of two variables, X and Y, we substitute X by X+Y:

Let’s open the first bracket (first binomic formula):

Now pull the “mean symbol” to each part of the expression:

Get rid of some brackets:

And we recognize some known expressions:

If X and Y are independent (uncorrelated) then Cov(X,Y) = 0. So we end up with:

Debrief

Well, we have dived in to the fact that the variance is additive. That is the variance of a sum (or difference) of k variables equals the sum of the variance of each variance (see expression above).

To come to this result we have used some simulation which has the advantage that it is less abstract. Then we flipped to some more rigorous proof (but which is less tangible).

What do we need all that stuff for? Basically a number of of statistical tests or procedures use this property such as the ANOVA or R squared. Maybe equally of more important, we here have an answer why the variance is used for descriptive statistics where some other measure appears more tangible. As the variance has this nice feature which we can use for later steps in data analytic projects, we compute it right at the beginning.

In dieser YACSDA (Yet-another-case-study-on-data-analysis) geht es um die beispielhafte Analyse nominaler Daten anhand des “klassischen” Falls zum Untergang der Titanic. Eine Frage, die sich hier aufdrängt, lautet: Kann (konnte) man sich vom Tod freikaufen, etwas polemisch formuliert. Oder neutraler: Hängt die Überlebensquote von der Klasse, in der derPassagiers reist, ab?

Diese Übung soll einige grundlegende Vorgehensweise der Datenanalyse verdeutlichen; Zielgruppe sind Einsteiger (mit Grundkenntnissen in R) in die Datenanalyse.

Daten laden

Zuerst laden wir die Daten. Es gibt mehrere Methoden, wie man Daten in R importieren kann. Eine einfache Möglichkeit ist, “Packages”, Pakete, zu nutzen. Einige Datensätze “wohnen” in R-Paketen. Man installiert also das entsprechende Paket, lädt das Paket und lädt dann drittens den Datensatz:

# install.packages("titanic")
library("titanic")
data(titanic_train)

Man beachte, dass ein Paket nur einmalig zu installieren ist (wie jede Software). Dann aber muss das Paket bei jedem Starten von R wieder von neuem gestartet werden. Außerdem ist es wichtig zu wissen, dass das Laden eines Pakets nicht automatisch die Datensätze aus dem Paket lädt. Man muss das oder die gewünschten Pakete selber (mit data(...)) laden. Und: Der Name eines Pakets (z.B. titanic) muss nicht identisch sein mit dem oder den Datensätzen des Pakets (z.B. titanic_train).

Erster Blick

Werfen wir einen ersten Blick in die Daten:

str(titanic_train)
## 'data.frame':	891 obs. of  12 variables:
##  $ PassengerId: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Survived   : int  0 1 1 1 0 0 0 0 1 1 ...
##  $ Pclass     : int  3 1 3 1 3 3 1 3 3 2 ...
##  $ Name       : chr  "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ...
##  $ Sex        : chr  "male" "female" "female" "female" ...
##  $ Age        : num  22 38 26 35 35 NA 54 2 27 14 ...
##  $ SibSp      : int  1 1 0 1 0 0 0 3 0 1 ...
##  $ Parch      : int  0 0 0 0 0 0 0 1 2 0 ...
##  $ Ticket     : chr  "A/5 21171" "PC 17599" "STON/O2. 3101282" "113803" ...
##  $ Fare       : num  7.25 71.28 7.92 53.1 8.05 ...
##  $ Cabin      : chr  "" "C85" "" "C123" ...
##  $ Embarked   : chr  "S" "C" "S" "S" ...

Ein ähnliches, etwas schön gegliedertes Ergebnis erreichen wir so:

# install.packages("dplyr", dependencies = TRUE) # ggf. vorher installieren
library(dplyr) 
glimpse(titanic_train)
## Observations: 891
## Variables: 12
## $ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,...
## $ Survived    <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0,...
## $ Pclass      <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3,...
## $ Name        <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bra...
## $ Sex         <chr> "male", "female", "female", "female", "male", "mal...
## $ Age         <dbl> 22, 38, 26, 35, 35, NA, 54, 2, 27, 14, 4, 58, 20, ...
## $ SibSp       <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4,...
## $ Parch       <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1,...
## $ Ticket      <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "1138...
## $ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, ...
## $ Cabin       <chr> "", "C85", "", "C123", "", "", "E46", "", "", "", ...
## $ Embarked    <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", ...

Welche Variablen sind interessant?

Von 12 Variablen des Datensatzes interessieren uns offenbar Pclass und Survived; Hilfe zum Datensatz kann man übrigens mit help(titanic_train) bekommen. Diese beiden Variablen sind kategorial (nicht-metrisch), wobei sie in der Tabelle mit Zahlen kodiert sind. Natürlich ändert die Art der Codierung (hier als Zahl) nichts am eigentlichen Skalenniveau. Genauso könnte man “Mann” mit 1 und “Frau” mit 2 kodieren; ein Mittelwert bliebe genauso (wenig) aussagekräftig. Zu beachten ist hier nur, dass sich manche R-Befehle verunsichern lassen, wenn nominale Variablen mit Zahlen kodiert sind. Daher ist es oft besser, nominale Variablen mit Text-Werten zu benennen (wie “survived” vs. “drowned” etc.). Wir kommen später auf diesen Punkt zurück.

Univariate Häufigkeiten

Bevor wir uns in kompliziertere Fragestellungen stürzen, halten wir fest: Wir untersuchen zwei nominale Variablen. Sprich: wir werden Häufigkeiten auszählen. Häufigkeiten (und relative Häufigkeiten, also Anteile oder Quoten) sind das, was uns hier beschäftigt.

Zählen wir zuerst die univariaten Häufigkeiten aus: Wie viele Passagiere gab es pro Klasse? Wie viele Passagiere gab es pro Wert von Survived (also die überlebten bzw. nicht überlebten)?

c1 <- count(titanic_train, Pclass)
c1
## # A tibble: 3 × 2
##   Pclass     n
##    <int> <int>
## 1      1   216
## 2      2   184
## 3      3   491

Aha. Zur besseren Anschaulichkeit können wir das auch plotten (ein Diagramm dazu malen). Ich empfehle den Befehl qplot, der mit konsistenter Syntax eine Vielzahl von Diagrammen ermöglicht.

# install.packages("ggplot2", dependencies = TRUE)
library(ggplot2)
qplot(x = Pclass, y = n, data = c1)

plot of chunk unnamed-chunk-5

Der Befehl qplot zeichnet automatisch Punkte, wenn auf beiden Achsen “Zahlen-Variablen” stehen (also Variablen, die keinen “Text”, sondern nur Zahlen beinhalten. In R sind das Variablen vom Typ int (integer), also Ganze Zahlen oder vom Typ num (numeric), also Reelle Zahlen).

c2 <- count(titanic_train, Survived)
c2
## # A tibble: 2 × 2
##   Survived     n
##      <int> <int>
## 1        0   549
## 2        1   342

Man beachte, dass der Befehl count stehts eine Tabelle (data.frame bzw. tibble) zurückliefert.

Bivariate Häufigkeiten

OK, gut. Jetzt wissen wir die Häufigkeiten pro Wert von Survived (dasselbe gilt für Pclass). Eigentlich interessiert uns aber die Frage, ob sich die relativen Häufigkeiten der Stufen von Pclass innerhalb der Stufen von Survived unterscheiden. Einfacher gesagt: Ist der Anteil der Überlebenden in der 1. Klasse größer als in der 3. Klasse?

Zählen wir zuerst die Häufigkeiten für alle Kombinationen von Survived und Pclass:

c3 <- count(titanic_train, Survived, Pclass)
c3
## Source: local data frame [6 x 3]
## Groups: Survived [?]
## 
##   Survived Pclass     n
##      <int>  <int> <int>
## 1        0      1    80
## 2        0      2    97
## 3        0      3   372
## 4        1      1   136
## 5        1      2    87
## 6        1      3   119

Da Pclass 3 Stufen hat (1., 2. und 3. Klasse) und innerhalb jeder dieser 3 Klassen es die Gruppe der Überlebenden und der Nicht-Überlebenden gibt, haben wir insgesamt 3*2=6 Gruppen.

Es ist hilfreich, sich diese Häufigkeiten wiederum zu plotten; wir nehmen den gleichen Befehl wie oben.

qplot(x = Pclass, y = n, data = c3)

plot of chunk unnamed-chunk-8

Hm, nicht so hilfreich. Schöner wäre, wenn wir (farblich) erkennen könnten, welcher Punkt für “Überlebt” und welcher Punkt für “Nicht-Überlebt” steht. Mit qplot geht das recht einfach: Wir sagen der Funktion qplot, dass die Farbe (color) der Punkte den Stufen von Survived zugeordnet werden sollen:

qplot(x = Pclass, y = n, color = Survived, data = c3)

plot of chunk unnamed-chunk-9

Viel besser. Was noch stört, ist, dass Survived als metrische Variable verstanden wird. Das Farbschema lässt Nuancen, feine Farbschattierungen, zu. Für nominale Variablen macht das keinen Sinn; es gibt da keine Zwischentöne. Tot ist tot, lebendig ist lebendig. Wir sollten daher der Funktion sagen, dass es sich um nominale Variablen handelt:

qplot(x = factor(Pclass), y = n, color = factor(Survived), data = c3)

plot of chunk unnamed-chunk-10

Viel besser. Jetzt noch ein bisschen Schnickschnack:

qplot(x = factor(Pclass), y = n, color = factor(Survived), data = c3) + 
  labs(x = "Klasse", 
       title = "Überleben auf der Titanic",
       colour = "Überlebt?")

plot of chunk unnamed-chunk-11

Signifikanztest

Manche Leute mögen Signifikanztests. Ich persönlich stehe ihnen kritisch gegenüber, da ein p-Wert eine Funktion der Stichprobengröße ist und außerdem zumeist missverstanden wird (er gibt nicht die Wahrscheinlichkeit der getesteten Hypothese an, was die Frage aufwirft, warum er mich dann interessieren sollte). Aber seisdrum, berechnen wir mal einen p-Wert. Es gibt mehrere statistische Tests, die sich hier potenziell anböten (was die Frage nach der Objektivität von statistischen Tests in ein ungünstiges Licht rückt). Nehmen wir den -Test.

chisq.test(titanic_train$Survived, titanic_train$Pclass)
## 
## 	Pearson's Chi-squared test
## 
## data:  titanic_train$Survived and titanic_train$Pclass
## X-squared = 102.89, df = 2, p-value < 2.2e-16

Der p-Wert ist kleiner als 5%, daher entscheiden wir uns gegen die H0 und für die H1: “Es gibt einen Zusammenhang von Überlebensrate und Passagierklasse”.

Effektstärke

Abgesehen von der Signifikanz, und interessanter, ist die Frage, wie sehr die Variablen zusammenhängen. Für Häufigkeitsanalysen mit 2*2-Feldern bietet sich das “Odds Ratio” (OR), das Chancenverhältnis an. Das Chancen-Verhältnis beantwortet die Frage: “Um welchen Faktor ist die Überlebenschance in der einen Klasse größer als in der anderen Klasse?”. Eine interessante Frage, als schauen wir es uns an.

Das OR ist nur definiert für 2*2-Häufigkeitstabellen, daher müssen wir die Anzahl der Passagierklassen von 3 auf 2 verringern. Nehmen wir nur 1. und 3. Klasse, um den vermuteten Effekt deutlich herauszuschälen:

t2 <- filter(titanic_train, Pclass != 2)  # "!=" heißt "nicht"

Alternativ (synonym) könnten wir auch schreiben:

t2 <- filter(titanic_train, Pclass == 1 | Pclass == 3)  # "|" heißt "oder"

Und dann zählen wir wieder die Häufigkeiten aus pro Gruppe:

c4 <- count(t2, Pclass)
c4
## # A tibble: 2 × 2
##   Pclass     n
##    <int> <int>
## 1      1   216
## 2      3   491

Schauen wir nochmal den p-Wert an, da wir jetzt ja mit einer veränderten Datentabelle operieren:

chisq.test(t2$Survived, t2$Pclass)
## 
## 	Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  t2$Survived and t2$Pclass
## X-squared = 95.893, df = 1, p-value < 2.2e-16

Ein Chi^2-Wert von ~96 bei einem n von 707.

Dann berechnen wir die Effektstärke (OR) mit dem Paket compute.es (muss ebenfalls installiert sein).

library(compute.es)
chies(chi.sq = 96, n = 707)
## Mean Differences ES: 
##  
##  d [ 95 %CI] = 0.79 [ 0.63 , 0.95 ] 
##   var(d) = 0.01 
##   p-value(d) = 0 
##   U3(d) = 78.59 % 
##   CLES(d) = 71.23 % 
##   Cliff's Delta = 0.42 
##  
##  g [ 95 %CI] = 0.79 [ 0.63 , 0.95 ] 
##   var(g) = 0.01 
##   p-value(g) = 0 
##   U3(g) = 78.56 % 
##   CLES(g) = 71.21 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = 0.37 [ 0.3 , 0.43 ] 
##   var(r) = 0 
##   p-value(r) = 0 
##  
##  z [ 95 %CI] = 0.39 [ 0.31 , 0.46 ] 
##   var(z) = 0 
##   p-value(z) = 0 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 4.21 [ 3.15 , 5.61 ] 
##   p-value(OR) = 0 
##  
##  Log OR [ 95 %CI] = 1.44 [ 1.15 , 1.73 ] 
##   var(lOR) = 0.02 
##   p-value(Log OR) = 0 
##  
##  Other: 
##  
##  NNT = 3.57 
##  Total N = 707

Die Chance zu überleben ist also in der 1. Klasse mehr als 4 mal so hoch wie in der 3. Klasse. Money buys you live…

Effektstärken visualieren

Zum Abschluss schauen wir uns die Stärke des Zusammenhangs noch einmal graphisch an. Wir berechnen dafür die relativen Häufigkeiten pro Gruppe (im Datensatz ohne 2. Klasse, der Einfachheit halber).

c5 <- count(t2, Pclass, Survived)
c5$prop <- c5$n / 707
c5
## Source: local data frame [4 x 4]
## Groups: Pclass [?]
## 
##   Pclass Survived     n      prop
##    <int>    <int> <int>     <dbl>
## 1      1        0    80 0.1131542
## 2      1        1   136 0.1923621
## 3      3        0   372 0.5261669
## 4      3        1   119 0.1683168

Genauer gesagt haben die Häufigkeiten pro Gruppe in Bezug auf die Gesamtzahl aller Passagiere berechnet; die vier Anteile addieren sich also zu 1 auf.

Das plotten wir wieder

qplot(x = factor(Pclass), y = prop, fill = factor(Survived), data = c5, geom = "col")

plot of chunk unnamed-chunk-19

Das geom = "col" heißt, dass als “geometrisches Objekt” dieses Mal keine Punkte, sondern Säulen (colons) verwendet werden sollen.

qplot(x = factor(Pclass), y = prop, fill = factor(Survived), data = c5, geom = "col")

plot of chunk unnamed-chunk-20

Ganz nett, aber die Häufigkeitsunterscheide von Survived zwischen den beiden Werten von Pclass stechen noch nicht so ins Auge. Wir müssen es anders darstellen.

Hier kommt der Punkt, wo wir von qplot auf seinen großen Bruder, ggplot wechseln sollten. qplot ist in Wirklichkeit nur eine vereinfachte Form von ggplot; die Einfachheit wird mit geringeren Möglichkeiten bezahlt. Satteln wir zum Schluss dieser Fallstudie also um:

ggplot(data = c5) +
  aes(x = factor(Pclass), y = n, fill = factor(Survived)) + 
  geom_col(position = "fill") +
  labs(x = "Passagierklasse", fill = "Überlebt?", caption = "Nur Passagiere, keine Besatzung")

plot of chunk unnamed-chunk-21

Jeden sehen wir die Häufigkeiten des Überlebens bedingt auf die Passagierklasse besser. Wir sehen auf den ersten Blick, dass sich die Überlebensraten deutlich unterscheiden: Im linken Balken überleben die meisten; im rechten Balken ertrinken die meisten.

Diese letzte Analyse zeigt deutlich die Kraft von (Daten-)Visualisierungen auf. Der zu untersuchende Effekt tritt hier am stärken zu Tage; außerdem ist die Analyse relativ einfach.

Fazit

In der Datenanalyse (mit R) kommt man mit wenigen Befehlen schon sehr weit; dplyr und ggplot2 zählen (zu Recht) zu den am häufigsten verwendeten Paketen. Beide sind flexibel, konsistent und spielen gerne miteinander. Die besten Einblicke haben wir aus deskriptiver bzw. explorativer Analyse (Diagramme) gewonnen. Signifikanztests oder komplizierte Modelle waren nicht zentral. In vielen Studien/Projekten der Datenanalyse gilt ähnliches: Daten umformen und verstehen bzw. “veranschaulichen” sind zentrale Punkte, die häufig viel Zeit und Wissen fordern. Bei der Analyse von nominalskalierten sind Häufigkeitsauswertungen ideal; ich hoffe, dieser Post konnte Ihnen weiterhelfen.