R-Pakete (packages) sind Erweiterungen für die statistische Programmiersprache R.
R-Pakete können Code, Daten und Dokumentation enthalten.

Sie können Pakete in der R Console mit dem Befehl install.packages installieren. Als Argument muss der Name des Packages in Anführungszeichen in der Klammer stehen.

Installieren Sie das Paket "palmerpenguins" mit folgendem Befehl:

install.packages("palmerpenguins")

Die Warnung: WARNING: Rtools is required to build R packages können Sie in diesem Fall ignorieren.

Achten Sie auf die Meldung: Paket ‘palmerpenguins’ erfolgreich ausgepackt und MD5 Summen abgeglichen diese zeigt an, dass alles passt.

Das Paket enthält gesammelte Echtdaten, die zum Üben verwendet werden können. Im Paket finden wir folgende Daten:

An 3 Pinguin-Arten (Adelie-Pinguine, Chinstrap-Pinguine, Gentoo-Pinguine) wurden an 344 Individuen verschiedene Körpermerkmale gemessen, z.B.

• Körpermasse (body mass) in Gramm (Spalte 'body_mass_g')
• Flossenlänge (flipper length) in Millimeter (Spalte 'flipper_length_mm')
• Schnabellänge (bill length) in Millimeter (Spalte 'bill_length_mm')

Details zu diesem Datensatz finden Sie unter: https://allisonhorst.github.io/palmerpenguins/

Öffnen Sie ein R Script Fenster:

Dieses können Sie dann gleich als .R-Dokument speichern, damit sie später auf den Code zugreifen, und Ihre Analysen wiederholen können:

Schreiben Sie in das Skript folgenden Befehl:

data('penguins', package='palmerpenguins')

und lassen Sie den Befehl mittels 'Run' (oder Tastenkombination Strg+Enter) ausführen:

Nun soll ein Histogramm mit den Werten der gemessenen Körpermassen (body_mass_g) der Pinguine gezeichnet werden:

hist(penguins$body_mass_g)

Führen Sie den Befehl wieder mittels 'Run' (oder Tastenkombination Strg+Enter) aus:

Beachten Sie, dass nur die Zeile ausgeführt wird, in der der Cursor steht!

Die Körpermasse aller gemessenen, erwachsenen Pinguine liegt also zwischen 2.5-6.5kg, wobei die meisten zwischen 3.5-4kg schwer sind. Die Verteilung ist leicht rechtsschief.

Nun soll auf der x-Achse die Flügellänge (flipper_length_mm) und auf der y-Achse die Körpermasse (body_mass_g) aufgetragen werden:

plot(penguins$flipper_length_mm, penguins$body_mass_g) 

In unserer Stichprobe besteht also ein positiver Zusammenhang zwischen Flügellänge und Körpermasse.

Lässt sich dieser Zusammenhang auch in der Grundgesamtheit wiederfinden? Um das zu beantworten, berechnen wir ein einfaches lineares Regressionsmodell. Dieses wird in der Variablen mod1 gespeichert, und mit der Funktion summary ausgegeben.

Das Tilde-Zeichen ~ geben Sie folgendermaßen ein:
Windows: AltGr plus Taste rechts neben dem 'ü'
Mac: siehe Mac Sonderzeichen

mod1 <- lm(body_mass_g ~ flipper_length_mm, data=penguins)
summary(mod1)

An der Ausgabe erkennen wir, dass sich knapp 76% der Gesamtvarianz in der Körpermasse der Pinguine aus unserer Stichprobe aus der Flügellänge erklären lassen. Das berechnete Regressionsmodell kann folgendermaßen aufgeschrieben werden:

$$ \text{body_mass_g} = -5780.8 + 49.7 \times \text{flipper_length_mm} $$, wobei sowohl der Achsenabstand als auch die Steigung der Grundgesamtheit bei einem Signifikanzniveau von $\alpha=0.05$ signifikant von 0 verschieden sind.

Wir können auch visuell die Annahmen für die Regression überprüfen:

par(mfrow=c(2,2))
plot(mod1)
par(mfrow=c(1,1))

Die Annahmen für die lineare Regression (Linearität und Normalverteilung der Fehler, Varianzgleicheit) werden also nicht verletzt, ebenso gibt es keine Ausreißer mit starkem Einfluss auf das Modell. Wir können also auf die Ergebnisse unseres Modells getrost vertrauen.

Nun wollen wir uns anschauen, ob denn der Zusammenhang (d.h. die Steigung der Regressionsgerade) zwischen Flügellänge und Körpermasse für alle 3 Spezies gleich ist.

Dafür plotten wir den Zusammenhang zuerst einmal in unterschiedlichen Farben für die drei Spezies:

plot(penguins$flipper_length_mm, penguins$body_mass_g,
     col=penguins$species)

Die Grafik lässt vermuten, dass das Verhältnis zwischen Flügellänge und Körpermasse nicht für alle Spezies gleich ist.

Wir überprüfen das in zwei weiteren Modellen:

mod2 <- lm(body_mass_g ~ flipper_length_mm + species, data=penguins)
summary(mod2)

Das resultierende Modell lässt sich also folgendermaßen aufschreiben: $$\text{body_mass_g} = -4031.5 + 40.7\times\text{flipper_length_mm}\\ - 206.5\times\text{species}_{\text{Chinstrap}} + 266.8\times\text{species}_{\text{Gentoo}}$$

Wo ist denn nun die Spezies 'Adelie' hinverschwunden? Tatsächlich werden Faktoren in R defaultmäßig so codiert, dass wir von einer Referenzkategorie ausgehen, mit der wir alle anderen vergleichen (um sich redundantes Rechnen zu ersparen). Defaultmäßig wird die Spezies als Referenz verwendet, die im Alphabet weiter vorne steht - R startet also mit Adelie. Die berechneten Werte für die anderen Spezies werden 'dazugeschaltet', wenn jeweils $\text{species}_{Chinstrap}=1$ oder $\text{species}_{Gentoo}=1$. Somit ist -4031.5 der Achsenabstand für Adelie, -4031.47-206.51=-4237.9 der Achsenabstand für Spezies Chinstrap, und -4031.477+266.81=-3764.667 der Achsenabstand für die Spezies Gentoo.

Visuell schaut das Modell folgendermaßen aus: Die Regressionslinien können Sie mit den Befehlen abline(coef=coefficients(mod2)[1:2]), abline(coef=c(sum(coefficients(mod2)[c(1,3)]), coefficients(mod2)[2]), col=2) und abline(coef=c(sum(coefficients(mod2)[c(1,4)]), coefficients(mod2)[2]), col=3) hinzufügen.

Das $R^2$ hat sich also durch das Berechnen eines separaten Achsenabstands pro Spezies auf 0.782 erhöht.

Können wir das Modell weiter verbessern, indem wir auch die Steigungen individuell berechnen?

mod3 <- lm(body_mass_g ~ flipper_length_mm * species, data=penguins)
summary(mod3)

Das neue Modell ist: $$\text{body_mass_g} = -2535.84 + 32.82\times\text{flipper_length_mm} \\ - 501.36 \times\text{species}_{\text{Chinstrap}} - 4251.44 \times\text{species}_{\text{Gentoo}} \\ + 1.74 \times \text{flipper_length_mm} \times \text{species}_{\text{Chinstrap}} \\ + 21.79 \times \text{flipper_length_mm} \times \text{species}_{\text{Gentoo}}$$

Nun sind analog zu oben der erste Achsenabstand und die erste Steigung die Werte für die Spezies Adelie, und Achsenabstand und Steigung für Chinstrap oder Gentoo lassen sich analog zu oben berechnen (z.B. die Steigung für Spezies Chinstrap ist nun 32.83+1.74=34.57).

Das vollständige Modell sieht folgendermaßen aus: An der Grafik spiegeln sich auch die Signifikanzen der Koeffizientenschätzer wider: Weder für den Achsenabstand, noch für die Steigung sind die Unterschiede zwischen Adelie und Chinstrap signifikant. Die Geraden sehen auch sehr ähnlich aus, und hätten wir Konfidenzbänder für unsere Regressionsgeraden gezeichnet, würden die beiden Geraden im Konfidenzband des jeweils anderen liegen. Die Unterschiede zwischen Adelie und Gentoo sind aber in beiden Fällen signifikant, und auch die berechente Gerade sieht vollkommen anders aus.

Der Anteil der erklärten Varianz erhöht sich um nur 0.007 im Vergleich zu mod2. Zahlt es sich also aus, diese zusätzlichen Variablen in das Modell aufzunehmen? Das testen wir erstmal mithilfe einer ANOVA aller 3 Modelle:

anova(mod1, mod2, mod3)

Sowohl der Anstieg um 3% im Anteil der erklärten Varianz zwischen mod1 und mod2, als auch der zusätzliche Anstieg um 0.7% zwischen mod2 und mod3 sind in der Grundgesamtheit signifikant von 0 verschieden, auch wenn letzterer in der Anwendung nur mehr sehr gering ist.

Auf die gleiche Weise könnten wir auch das Hinzufügen anderer erklärender Variablen testen. Um uns hier aber Tipp-Arbeit zu ersparen, nutzen wir die Selektionsfuntion step(), die anhand des Akaike Information Criterion (AIC) Variablen abwechselnd hinzufügt und wegnimmt, und die Veränderung der Erklärkraft des Modells testet:

mod4 <- lm(body_mass_g~., data=penguins)
mod5 <- step(mod4, direction='both')
summary(mod5)

Dabei haben wir zuerst mod4 als Startmodell erstellt. Dabei heißt '~.', dass 'body_mass_g' aus allen anderen Variablen im Datensatz erklärt wird. (Ebenso könnten wir mit '~.^2' verlangen, dass, so wie in mod3 im Vergleich zu mod2, statt einem + ein * zwischen die Variablen gesetzt wird, also auch alle möglichen Interaktionen zwischen Variablen getestet werden. Das ersparen wir uns aber aus Zeitgründen.) Danach werden die Variablen des großen Modells mod4 in der step()-Funktion selektiert, und das resultierende Modell mod5 wird angezeigt. (Auch der Befehl step() selbst produziert Output an der Konsole, da dieser aber sehr lang ist, sparen wir uns den hier.)

Laut AIC ist also das beste Modell jenes, bei dem die Körpermasse aus Spezies, Schnabellänge, Schnabeltiefe, Flügellänge, Geschlecht, und Jahr der Messung erklärt wird. Das $R^2$ des Modells ist 0.876, wir können also knapp 88% der Varianz in der Körpermasse der Pinguine aus unserem Modell erklären.