Auf dieser Seite findest du Übungen zu Multipler Linearer Regression in R.

Der gesamte Code, den du hier findest, kann von dir 1:1 in R umgesetzt werden. Die größte Lernwirkung hat diese Einführung, wenn du möglichst alles einmal selbst in R durchführst.

01 Vorbereitungen

Bevor wir loslegen können, müssen wir die Packages laden, die wir für unser Vorhaben benötigen. Zunächst ist das nur das SfL Package. Falls du dieses noch nicht geladen hast, hole dies nun nach:

library(SfL)

Nun laden wir die Datensätze, die wir nutzen werden:

data(data_a)
data(data_c)
data(data_s)
data(data_t)
data(data_v)

02 Einführung

Multiple Lineare Regression wird ebenso wie Einfache Lineare Regression wird in R mit dem lm() Befehl umgesetzt. Dieser wird genutzt, um zu definieren, welche Variablen voneinander abhängig modelliert werden sollen.

Multiple Lineare Regression ist eine Erweiterung der allgemeinen Formel für Einfache Lineare Regression:

\[Y=\beta_1+\beta_2X+\epsilon\]

Hierbei beschreibt \(Y\) die Variable, welche wir modellieren, d.h. vorhersagen, möchten. \(X\) ist die Variable mit welcher wir \(Y\) vorhersagen möchten.

Hier nun die allgemeine Formel für Multiple Lineare Regression:

\[Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+...+\beta_{i}X_{i}+\epsilon\]

Wie zuvor ist auch hier \(Y\) die Variable, welche wir modellieren möchten. Anders als zuvor haben wir nun mehrere - multiple - \(X\) Variablen im Modell. All diese Variablen werden genutzt, um \(Y\) vorherzusagen. Jede dieser \(X\) Variablen verfügt über einen eigenen eigenen Slope: \(\beta_{i}\). Der Intercept des Modells ist durch \(\beta_{0}\) angegeben. \(\epsilon\) ist wie zuvor auch der Fehlerterm des Modells.

03 Modelling

Als Beispiel erstellen wir nun ein Einfaches Lineares Modell für die Dauer von /s/. Vorhergesagt werden soll die Dauer durch die Sprechgeschwindigkeit und die log-transformierte Dauer der Base. Wichtig ist auch bei Regressionsmodellen, dass die abhängige Variable möglichst normalverteilt ist. Daher überprüfen wir dies im nächsten Schritt.

Distribution Check

Wie wir bereits in vorherigen Sessions gesehen haben, ist sDur nicht normalverteilt:

shapiro.test(data_s$sDur)

    Shapiro-Wilk normality test

data:  data_s$sDur
W = 0.93622, p-value = 2.761e-06

Wir nutzen daher eine log-transformierte Version der Variable als abhängige Variable im weiteren Vorgehen:

data_s$sDurLog <- log(data_s$sDur)

Modell-Erstellung

Nun erstellen wir das Modell:

lm_02 <- lm(sDurLog ~ speakingRate + baseDurLog, data = data_s)

Regressionsmodelle sollten immer als Objekt gespeichert werden. Daher nennen wir dieses erste Modell lm_02. So ist sichergestellt, dass wir das Modell auch mit weiteren Funktionen nutzen können.

Koeffizienten

Nun schauen wir uns zunächst an, was R ausgibt, wenn wir das erstellte Modell abfragen:

lm_02

Call:
lm(formula = sDurLog ~ speakingRate + baseDurLog, data = data_s)

Coefficients:
 (Intercept)  speakingRate    baseDurLog  
    -1.03540      -0.03138       0.82203

Hier finden wir erneut die Formel zu unserem Modell sowie das genutze Datenset. Außerdem erhalten wir die Coefficients. Doch was beudeuten diese Angaben?

  • (Intercept): der Y-Achsenabschnitt unserer Regressionsgeraden, d.h. die durchschnittliche vorhergesagte Dauer eines /s/ bei niedrigster Sprechgeschwindigkeit.
  • speakingRate: Entspricht \(X_{1}\) in der generellen Formel; dies ist der Regressionskoeffizient der Variable.
  • baseDurLog: Entspricht \(X_{2}\) in der generellen Formel; dies ist der Regressionskoeffizient der Variable.

Signifikanz

Interessant nach dem Erstellen eines Modells ist natürlich die Frage, ob \(Y\) wirklich von \(X\) vorhergesagt werden kann, also ob im vorliegenden Fall die Dauer von /s/ durch die Sprechgeschwindigkeit und/oder die Base Duration modelliert wird. Dies finden wir mit der anova() Funktion heraus:

anova(lm_02)
Analysis of Variance Table

Response: sDurLog
              Df  Sum Sq Mean Sq F value    Pr(>F)    
speakingRate   1  1.8858  1.8858  15.921 0.0001038 ***
baseDurLog     1  5.6391  5.6391  47.609 1.442e-10 ***
Residuals    147 17.4115  0.1184                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hier erhalten wir einige Werte zu den Bestandteilen des Modells. Für uns interessant ist der p-Wert (Pr(>F)): liegt dieser unter \(0.05\), so können wir davon ausgehen, dass \(X\) tatsächlich \(Y\) beeinflusst, d.h. modelliert.

Finding the Best Model

Nutzen wir Simple Lineare Regression, so ist jedes Modell ein finales Modell. Da wir nur eine Variable als vorhersagende Variable nutzen können, gibt es auch keine Auswahl zu treffen bzgl. der Variablenwahl.

Arbeiten wir allerdings mit Multipler Linearer Regression, so können wir theoretisch unendlich viele Variablen in ein Modell aufnehmen. Allerdings stellt sich oftmals heraus, dass das Behalten von Variablen ohne signifikante Effekte ein Modell insgesamt verschlechtert. Daher stellen sich die Fragen: Welche Variablen haben einen Effekt? Welche Variablen schaffen es in mein finales Modell?

Um diese Fragen zu beantworten, nehmen wir folgende zwei Punkte an:

  1. Wir haben erneut vor, sDurLog zu modellieren
  2. Wir haben Daten zu folgenden Dingen:
    speech rate, type of /s/, base duration, following pause, and biological sex

Mit diesen zwei Punkten können wir ein s.g. volles Modell erstellen:

lm_full <- lm(sDurLog ~ speakingRate + typeOfS + baseDurLog + pauseBin + bsex, data = data_s)

Schauen wir uns das volle Modell genauer an:

anova(lm_full)
Analysis of Variance Table

Response: sDurLog
              Df  Sum Sq Mean Sq F value    Pr(>F)    
speakingRate   1  1.8858  1.8858 21.7477 7.068e-06 ***
typeOfS        2  3.8115  1.9058 21.9779 4.755e-09 ***
baseDurLog     1  4.6265  4.6265 53.3544 1.790e-11 ***
pauseBin       1  1.9002  1.9002 21.9138 6.557e-06 ***
bsex           1  0.3126  0.3126  3.6049   0.05962 .  
Residuals    143 12.3998  0.0867                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anscheinend sind alle Variablen außer bsex signifikante Prediktoren für sDurLog.

Nun zurück zum eigentlichen Problem: Wie finden wir das beste Modell? Eine Möglichkeit wäre es, manuell einzelne Variablen zu entfernen und sich dabei an aufgestellten Kriterien zu orientieren. Diese Möglichkeit ist allerdings zeitaufwendig und fehleranfällig. Zum Glück stellt uns R eine viel bequemere Möglichkeit zur Verfügung:

step(lm_full)
Start:  AIC=-359.94
sDurLog ~ speakingRate + typeOfS + baseDurLog + pauseBin + bsex

               Df Sum of Sq    RSS     AIC
- speakingRate  1   0.14962 12.549 -360.14
<none>                      12.400 -359.94
- bsex          1   0.31259 12.712 -358.21
- pauseBin      1   1.75782 14.158 -342.06
- typeOfS       2   2.81457 15.214 -333.26
- baseDurLog    1   2.67066 15.070 -332.68

Step:  AIC=-360.14
sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex

             Df Sum of Sq    RSS     AIC
<none>                    12.549 -360.14
- bsex        1    0.3856 12.935 -357.60
- pauseBin    1    1.7686 14.318 -342.37
- typeOfS     2    2.7324 15.282 -334.60
- baseDurLog  1    3.7772 16.327 -322.68

Call:
lm(formula = sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex, 
    data = data_s)

Coefficients:
  (Intercept)      typeOfSnm      typeOfSpl     baseDurLog  pauseBinpause          bsexm  
      -1.5954         0.3266         0.1015         0.6665         0.2390         0.1148

Die erste Reihe des Outputs zeigt den AIC-Wert des vollen Modells an: \(-359.94\). Mit dem AIC-Wert kann man die Güte eines Modells angeben: Je kleiner der Wert, desto besser das Modell.

Darunter finden wir die Formel des vollen Modells: sDurLog ~ speakingRate + typeOfS + baseDurLog + pauseBin + bsex.

Es folgt eine Tabelle, die alle Variablen des vollen Modells beinhaltet. Für jede Variable erstellt die step() Funktion ein Modell ohne eben diese Variable. Beispiel: Ein Modell mit allen Variablen außer bsex hat einen AIC-Wert von \(-358.21\).

Anhand dieses Vorehens findet die step() Funktion schließlich das Modell mit der besten Variablenkombination, d.h. das Modell mit dem niedrigsten AIC-Wert. Die Formel dieses Modells ist unter Call angegeben: sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex.

Zusätzlich wird uns der AIC-Wert des gefundenen Modells mitgeteilt; so können wir diesen direkt mit dem Ausgangswert vergleichen.

Zuletzt werden die Koeffizienten des neuen Modells gegeben.


Mit der gefundenen Formel können wir nun das finale Modell erstellen:

lm_final <- lm(formula = sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex, data = data_s)

Summary

Wie auch bei Simpler Linearer Regression kann auch hier die summary() Funktion genutzt werden:

summary(lm_final)

Call:
lm(formula = sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex, 
    data = data_s)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.77731 -0.17999  0.01121  0.20711  0.79500 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.59544    0.13904 -11.475  < 2e-16 ***
typeOfSnm      0.32658    0.05984   5.458 2.05e-07 ***
typeOfSpl      0.10147    0.05997   1.692   0.0928 .  
baseDurLog     0.66652    0.10124   6.583 8.00e-10 ***
pauseBinpause  0.23899    0.05305   4.505 1.36e-05 ***
bsexm          0.11482    0.05459   2.103   0.0372 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2952 on 144 degrees of freedom
Multiple R-squared:  0.4967,    Adjusted R-squared:  0.4793 
F-statistic: 28.43 on 5 and 144 DF,  p-value: < 2.2e-16

Tukey Contrasts

Tukey Contrasts werden genutzt um die Level von kategorischen Prediktoren zu untersuchen, d.h. um zu schauen, ob diese signifikant verschieden zueinander sind.

library(SfL)

tukey(model = lm_final, predictor = typeOfS)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex, 
    data = data_s)

Linear Hypotheses:
             Estimate Std. Error t value Pr(>|t|)    
nm - is == 0  0.32658    0.05984   5.458  < 1e-04 ***
pl - is == 0  0.10147    0.05997   1.692 0.211759    
pl - nm == 0 -0.22512    0.05945  -3.787 0.000662 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
tukey(model = lm_final, predictor = pauseBin)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex, 
    data = data_s)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)    
pause - no_pause == 0  0.23899    0.05305   4.505 1.36e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
tukey(model = lm_final, predictor = bsex)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = sDurLog ~ typeOfS + baseDurLog + pauseBin + bsex, 
    data = data_s)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)  
m - f == 0  0.11482    0.05459   2.103   0.0372 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

04 Diagnostik

Wie wir bereits im Theorieteil gelernt haben, sind lineare Modelle nicht verlässlich, wenn ihre Annahmen (assumptions) verletzt sind. Daher checken wir diese in den folgenden Abschnitten.

Linearity

library(performance)

plot(check_heteroscedasticity(lm_final))
OK: Error variance appears to be homoscedastic (p = 0.294).

The line should be horizontal and flat.

Homoscedasticity

library(performance)

plot(check_heteroscedasticity(lm_final))
OK: Error variance appears to be homoscedastic (p = 0.294).

Data points should be spread equally around the line, with no obvious patterns visible.

Normality

library(performance)

plot(check_normality(lm_final))
OK: residuals appear as normally distributed (p = 0.374).

The distribution of a linear model’s residuals should be normally distributed as indicated by the greenish/blueish curve.

Independence

Independence cannot be checked visually. It is an assumption that you can test by examining the study design.

05 Aufgaben

Erstelle nun Einfache Lineare Modelle für die folgenden \(X\)- und \(Y\)-Variablen. Sei dir bei jedem Modell darüber im Klaren, ob die vorhersagende Variable kontinuierlich oder kategorisch ist.

Aus data_s:

  • sDurLog ~ bsex + age + location
  • sDurLog ~ typeOfS + pauseBin
  • sDurLog ~ baseDurLog + speakingRate
  • sDurLog ~ biphoneProbSumBin + folType + preC
  • sDurLog ~ bsex + age + location + typeOfS + pauseBin + baseDurLog + speakingRate + biphoneProbSumBin + folType + preC


Aus data_a:

  • height ~ bsex + hair
  • height ~ hair + eyes
  • height ~ bsex + hair + hair + eyes


Aus data_v:

duration ~ vowel + structure + rate