Das allgemeine lineare Modell (Teil 1)

Übung Forschungsmethodik II | SS 2024

Autor
TU Chemnitz

C. Maiwald

Professur Forschungsmethoden und Analyseverfahren in der Biomechanik

Datum

11. April 2024

1 Ein lineares Modell

Lineare Modelle stellen einen Zusammenhang zwischen einer unabhängigen Variable (UV) und einer abhängigen Variable (AV) her. Lineare Modelle bzw. Zusammenhänge zweier Variablen lassen sich in der Form \[ AV = b_0 + b_1 \cdot UV\] beschreiben. Ein Mehr in der UV resultiert je nach Ausprägung des Modellparameters \(b_1\) in einem Mehr (wenn \(b_1\) negativ ist: Weniger) in der Vorhersage der AV. Damit entstehen bei konstanten Inkrementen der UV stets konstante/proportionale Veränderungen der AV. Im Gegensatz dazu wäre bei einem nicht-linearen Modell die Veränderung der AV über den Wertebereich der UV nicht konstant, sondern variierte je nach Ausprägung bzw. Ausgangszustand der UV.

Wir werden sehen, dass diese Art des linearen Modells (die allgemeinste Form, das allgemeine lineare Modell) die Grundlage für zahlreiche bisher als unterschiedlich angesehene Testverfahren bildet und damit eine sehr weitreichende Bedeutung entfalten kann.

2 Der t-Test im ALM

Im ersten Beispiel entwickeln wir eine Alternative zum t-Test mittels ALM. Wir nutzen ein lineares Modell, um den Zusammenhang zwischen der dichotomen Gruppierungsvariable Bluthochdruck (UV) HYPTEN und der metrischen AV Bodymass-Index (BMI) zu analysieren. Dies ist – wie wir noch sehen werden – im Ergebnis identisch zu einer Durchführung eines t-Tests, bei dem wir prüfen, inwiefern sich die Mittelwerte des BMI zwischen den beiden Gruppen in HYPTEN unterscheiden.

2.1 Datenvorbereitung

Wir zeigen die Ausgangslage anhand des Datensatzes physicalActivity.RData (siehe OPAL). Wir bilden zunächst den Subdatensatz:

load("physicalActivity.RData")
D <- subset(physicalActivity, select=c(HYPTEN, BMI))

2.2 Anpassung des linearen Modells auf die Daten

Der Code zur Durchführung der linearen Modellierung ist in R denkbar einfach gehalten. Die Funktion lm() stellt das linear model zur Verfügung, welches wir nur noch mit der entsprechenden Konstellation der AV und UV in Formelschreibweise füllen müssen:

LM.BMI <- lm(BMI ~ HYPTEN, data=D)

Damit sind alle Berechnungsschritte der linearen Modellierung erledigt, inklusive der Anpassung der Modellparameter \(b_0\) und \(b_1\), so dass die Residuale (d. h. die Abweichungen der Datenpunkte von den Vorhersagen der Gerade) minimiert und in Mittel gleich Null sind.

Die Modellausgabe erzeugen wir mit der Funktion summary():

summary(LM.BMI)

Call:
lm(formula = BMI ~ HYPTEN, data = D)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.4719 -2.8705 -0.0362  2.8281 10.8995 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  26.6005     0.6522  40.787   <2e-16 ***
HYPTENYes     2.2713     0.9646   2.355   0.0214 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.02 on 68 degrees of freedom
Multiple R-squared:  0.07539,   Adjusted R-squared:  0.0618 
F-statistic: 5.545 on 1 and 68 DF,  p-value: 0.02143

Wir erkennen in der Ergebnisausgabe mehrere Unterabschnitte:

  • Call: nochmalige Darstellung der AV und UV in Formelschreibweise, zur Überprüfung und richtigen Zuordnung der Ausgabe.
  • Residuals: five-number-summary der Residuale
  • Coefficients: der für uns wichtige Teil, welcher die Koeffizienten \(b_0\) und \(b_1\) enthält (Spalte Estimate), sowie weitere Angaben zum Signifikanztest der Koeffizienten
  • Residual standard error: Der “im Mittel” (siehe VL-Skript) zu erwartende Betrag der Residuale. Je kleiner, desto besser gelingen die Vorhersagen mit Hilfe des Modells.
  • Multiple R-squared und Adjusted R-squared: Angaben zur Varianzaufklärung des Modells. Je näher die Werte an 1 liegen, desto besser funktioniert das Modell für die Vorhersage der y-Werte. Bei \(R^2=1\) wäre eine perfekte Vorhersage gegeben, d. h. 100% Varianzaufklärung durch das lineare Modell. Dazu müssten aber alle Datenpunkte auf der Geraden liegen und die Residuale gleich Null sein. Kommt in der Realität daher nie vor, man erreicht die 1 nur näherungsweise in besonders geeigneten Variablenkonstellationen.
  • F-statistic: die Signifikanzprüfung des gesamten Modells. Je höher der F-Wert, desto kleiner der p-Wert (vgl. Signifikanztest einer ANOVA: F-Test).

2.3 Vergleich zum t-Test

Wie oben angerissen könnten wir diese Analyse auch mit Hilfe eines t-Tests durchführen. Dazu prüfen wir den Mittelwertsunterschied des BMI zwischen den Gruppen:

T.BMI <- t.test(BMI ~ HYPTEN, data=D, var.equal=TRUE)
print(T.BMI)

    Two Sample t-test

data:  BMI by HYPTEN
t = -2.3547, df = 68, p-value = 0.02143
alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
95 percent confidence interval:
 -4.1961506 -0.3465468
sample estimates:
 mean in group No mean in group Yes 
         26.60053          28.87187 

Wir sehen, dass wir für diesen t-Test eine vom Betrag her identische Ausprägung des t-Werts, der Freiheitsgrade und des p-Werts erhalten. In der Liste der Ausgabe zu T.BMI versteckt sich auch der Standardfehler, der mit stderr = 0.9646 angegeben wird. Die im t-Test ermittelte Differenz der Gruppenmittelwerte versteckt sich als Differenz der beiden ausgegebenen Gruppenmittel in der Variable T.BMI$estimate und beläuft sich auch hier auf den Wert diff(T.BMI$estimate) = 2.2713.

Damit haben wir aus zwei vordergründig völlig unterschiedlichen Herangehensweisen exakt dieselben Ergebnisse erhalten:

  • Der Koeffizient \(b_0\) beschreibt den Vorhersagewert (d. h. Mittelwert) für den BMI in der Gruppe mit normalem Blutdruck (Referenzgruppe, Faktorlevels in HYPTEN beachten! ) und damit den Achsenabschnitt (die Konstante, engl. intercept) der Modellrepräsentation einer Gerade.
  • Der Koeffizient \(b_1\) beschreibt die Differenz der beiden Gruppenmittelwerte, im Falle einer Modellrepräsentation als Gerade gibt er die Steigung (engl. slope) an.
  • Die Summe der beiden Koeffizienten \(b_0 + b_1\) beschreibt den Vorhersagewert (d. h. Mittelwert) der Bluthochdruckgruppe.
  • Der t-Test des Mittelwertsunterschieds zwischen den Gruppen ist folglich identisch zum t-Test des Koeffizienten \(b_1\): beide erzeugen folglich auch dieselben Prüfstatistiken.
Vorsicht

Diese Interpretation setzt die übliche Codierung der Gruppen in HYPTEN mit den Werten 0 und 1 voraus.

2.4 Grafische Modellrepräsentation

Modellrepräsentation als Gerade bei nicht-metrischer x-Achse

Die Darstellung einer Gerade erfordert normalerweise zwei metrische Achsen im Schaubild. Bei der Darstellung einer Gerade als Modellrepräsentation des linearen Modells mit einer nicht-metrischen Gruppierungsvariable auf der x-Achse führt dies in R konsequenterweise zu Problemen. Wir zeigen aus illustrativen Zwecken hier, wie man diese Probleme umgehen kann. Dies dient lediglich dazu, eine bildliche Vorstellung über die identischen Aussagen von t-Test und ALM zu schaffen und ist natürlich inhaltlich höchst fragwürdig, da diese Gerade nur an den Punkten der beiden Gruppenmittelwerte sinnvoll interpretiert werden kann. Wir raten dringend davon ab, solch eine Darstellung in wissenschaftlichen Abhandlungen tatsächlich anzuwenden.

Die folgenden Prozeduren der grafischen Repräsentation eines linearen Modells bei nicht-metrischer x-Achse unterscheiden sich zwischen base R und ggplot2, wobei jeweils andere Probleme ursächlich und dementsprechend zu korrigieren sind.

2.4.1 base R

Um eine Vorstellung der Anpassung eines ALM auf die vorliegenden Daten zu bekommen, können wir uns die resultierende Gerade als zusätzliche Linie in eine bestehende Grafik plotten:

plot(BMI~HYPTEN, data=D)
abline(LM.BMI, col="red")

Vorsicht

Es ist offensichtlich, dass in der grafischen Darstellung eines base R Boxplots die Codierung der x-Achse nicht zur eingezeichneten Regressionsgerade passt. Die rote Gerade müsste eigentlich durch den Mittelwert der Gruppe No verlaufen – das tut sie erkennbar nicht. Wo liegt das Problem?

Wir sehen uns die Lage des Werts Null auf der x-Achse der Grafik an, indem wir den Wertebereich auf dieser Achse erweitern und eine gestrichelte vertikale Linie an der Stelle x=0 einfügen:

plot(BMI ~ HYPTEN, data=D, xlim=c(-1,3))
abline(LM.BMI, col="red")
abline(v=0, lty=2, col="blue")

Man erkennt, dass R die Position der Boxen bei der Erzeugung eines Boxplot (in base R) nicht automatisch mit 0 beginnt, sondern die Boxen beginnend beim Wert 1 platziert. Wie korrigiert man das?

Standardmäßig werden die Boxen in base-R Boxplots laut Hilfefunktion an den Positionen “1:n where n is the number of boxes.” platziert. Das Argument at=c() in der Funktion boxplot() kann diese Standardposition korrigieren, wir müssen also die Positionen der Boxen händisch verändern, wenn wir sie im Kontext eines ALM zur Position der Geraden (die standardmäßig auf die Werte 0 und 1 codiert sind) anpassen wollen:

boxplot(BMI ~ HYPTEN, data=D, xlim=c(-1,3), at=c(0,1))
abline(LM.BMI, col="red")
abline(v=0, lty=2, col="blue")

Vorsicht

Um die Position von base R Boxplots mit jenen der Modellausgabe eines ALM in Übereinstimmung zu bringen, ist die händische Korrektur der Boxpositionen erforderlich. Dies geschieht mit dem Funktionsargument at=c(), welches einen Vektor der gewünschten Positionen enthalten muss. Andernfalls beginnt die Reihenfolge der Boxen an der Position 1, und stimmt damit nicht mit der Codierung der Geraden überein, welche die erste Box auf Position 0 der x-Achse erwartet.

2.4.2 ggplot2

P <- ggplot(D, aes(x=HYPTEN, y=BMI))+
  geom_boxplot()+
  geom_smooth(method="lm")

P
`geom_smooth()` using formula = 'y ~ x'

In ggplot2 geschieht an dieser Stelle nichts, wir müssen zunächst die Ursache für die nicht eingezeichnete Gerade suchen.

Die Variable HYPTEN ist eine Faktorvariable mit zwei Levels. Aus der Hilfe zu geom_smooth() entnehmen wir, dass diese Visualisierung nur für numerische Variablen funktioniert. Owohl R Faktorvariablen im Hintergrund als ganzzahlige numerische Variablen repräsentiert, wird diese Codierung von ggplot2 offenbar nicht genutzt. Wir müssen sie also händisch erstellen.

D$HYPTEN.num <- as.integer(D$HYPTEN)

P <- ggplot(D, aes(x=HYPTEN, y=BMI))+
  geom_boxplot()+
  geom_smooth(mapping=aes(x=HYPTEN.num, y=BMI), method="lm", se=FALSE)

P
`geom_smooth()` using formula = 'y ~ x'

Einfacher mit ggplot2

Der Umstand, dass diese Methode ohne Korrektur der Box-Positionen auskommt, liegt daran, dass sowohl die Faktorlevels in R als auch die Boxpositionen in ggplot2 beim Wert 1 beginnen. Für eine korrekte Positionierung der Regressionsgerade muss also nur sichergestellt werden, dass eine numerische Version der Gruppierungsvariable (x-Achse) bereitsteht.

2.5 Grafische Modellprüfung

Einige der Anwendungsvoraussetzungen lassen sich im ALM erst nach der Anpassung des Modells auf die Daten beantworten. Insbesondere betrifft dies die Annahmen zur Verteilung der Residuale: sind diese zufällig und in etwa normal verteilt?

Wir prüfen diese Dinge effektiv mittels der plot() Funktion, welche uns automatisch die Diagnostic plots des linearen Modells liefert. Insgesamt handelt es sich um vier Plots, die für uns wichtigen Plots sind die ersten beiden:

par(mfrow=c(1,2)) # Erzeugen zweier Grafiken nebeneinander (Anordnung: Eine Zeile, zwei Spalten)
plot(LM.BMI, c(1:2))

Wir sehen am ersten Plot (Residual Plot), dass die Residuale (in etwa) gleichermaßen streuen, und kein erkennbarer Trend vorliegt (rote Linie annähernd horizontal).

Im Q-Q Plot erkennen wir eine annähernde Normalverteilung der Residuale anhand der Lage vieler Datenpunkte auf der gestrichelten Diagonale.

Damit sind die wichtigsten Kriterien zur Anwendung des ALM erfüllt, und wir können von belastbaren Ergebnissen unserer Modellierung ausgehen.

2.6 Ergebniszusammenfassung

Wir können im Ergebnis feststellen, dass wir einen statistisch signifikanten Zusammenhang zwischen dem Vorliegen eines Bluthochdrucks und BMI ermittelt haben. Das ist gleichbedeutend mit einem statistisch signifikanten Mittelwertsunterschied (Differenz: ca. 2.27 \(kg/m^2\)) im BMI zwischen beiden Gruppen. Durch die Reihenfolge der Faktorlevels erkennen wir ebenfalls die Effektrichtung: der Koeffizient \(b_1\) ist positiv, daher ist der BMI-Mittelwert der Bluthochdruck-Gruppe größer als jener der Normal-Gruppe.

Zwischenfazit: t-Test vs. ALM

Das ALM enthält sämtliche Informationen, die ein t-Test bereitstellen kann. Es bietet darüber hinaus eine größere Flexibilität in der Anwendung. Wir können den t-Test (für unabhängige Stichproben) daher durch die Anwendung des ALM ersetzen.

3 Die einfache lineare Regression im ALM

Das, was wir bisher unter dem Begriff der einfachen linearen Regression behandelt haben, ist letztendlich ebenfalls nur ein Spezialfall des ALM. Die Verwandschaft ist offensichtlich, denn die Gleichung des linearen Modells im Abschnitt oben unterscheidet sich letztendlich nicht von einer Regressionsgleichung der einfachen linearen Regression. Auch die Modellrepräsentation in Form einer Gerade entspricht dem, was wir als Regressionsgerade kennen.

3.1 Neues Modell mit metrischer UV (Prädiktor): WAIST.CIR

Wir sehen uns den Zusammenhang des BMI mit der UV WAIST.CIR an, dem Taillenumfang (in cm). Ein erhöhter Taillenumfang sollte mit einem höheren BMI einhergehen. Allerdings spielt auch die Körpergröße eine gewisse Rolle, denn größere Menschen haben per se einen höheren Taillenumfang, möglicherweise eben auch unabhängig vom BMI. Wir sollten daher den Taillenumfang relativ zur Körpergröße (in Prozent) ausdrücken. Je größer dieser prozentuale Wert wird, desto eher würden wir einen höheren BMI erwarten.

Wir passen ein neues Modell auf die neue Variablenkonstellation an. Im Ergebnis erhalten wir:

D$WAIST <- physicalActivity$WAIST.CIR/physicalActivity$HEIGHT * 100
LM.BMI2 <- lm(BMI ~ WAIST, data=D)
summary(LM.BMI2)

Call:
lm(formula = BMI ~ WAIST, data = D)

Residuals:
   Min     1Q Median     3Q    Max 
-4.465 -1.085  0.159  1.253  2.940 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.99717    1.57381  -1.269    0.209    
WAIST        0.52121    0.02746  18.983   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.666 on 68 degrees of freedom
Multiple R-squared:  0.8413,    Adjusted R-squared:  0.8389 
F-statistic: 360.4 on 1 and 68 DF,  p-value: < 2.2e-16

Da in R auch einfache lineare Regressionen mit derselben Funktion lm() berechnet werden, müssen wir keine Vergleiche anstellen: das ALM stellt uns direkt die Ausgabe einer einfachen linearen Regression bereit, sofern wir die Variablenkonstellation entsprechend wählen (AV: metrisch, eine UV: metrisch).

Die Regressionsgleichung (bzw. die Gleichung des ALM) mit den Koeffizienten \(b_0\) und \(b_1\) entnehmen wir der Ausgabe:

\[\hat{y}_i = -1.99717 + 0.52121 \cdot x_i \] ## Grafische Repräsentation des Modells Wir erzeugen uns eine geeignete Grafik, in die wir die Ausgaben des Modells integrieren. Dazu wählen wir die Darstellung als Streudiagramm.

3.1.1 base R

plot(BMI ~ WAIST, data=D)
abline(LM.BMI2, col="red")

Die Ausgabe des Modells (Geradengleichung) lässt sich in base R einfach durch die Funktion abline() in die Grafik einbauen. In diesem Fall (metrische Variable auf x-Achse) stimmt die Position der Geraden auch auf Anhieb.

3.1.2 ggplot2

P2 <- ggplot(D, aes(x=WAIST, y=BMI))+
  geom_point()+
  geom_smooth(method="lm", se=FALSE)

P2
`geom_smooth()` using formula = 'y ~ x'

Auch in ggplot2 erzeugt die Anweisung geom_smooth(method="lm") automatisch die korrekt positionierte Regressionsgerade in der Grafik.

3.1.3 Exkurs: Ergänzung der Grafik um Referenzbereiche

Die üblichen Referenzbereiche des BMI lassen sich zu illustrativen Zwecken recht einfach durch ribbons in der Grafik ergänzen:

P3 <- P2 + 
  theme_minimal()+
  geom_ribbon(aes(ymin=18, ymax=25), fill="green", alpha=0.2)+
  geom_ribbon(aes(ymin=25, ymax=30), fill="orange", alpha=0.2)+
  geom_ribbon(aes(ymin=30, ymax=35), fill="red", alpha=0.2)+
  geom_ribbon(aes(ymin=35, ymax=40), fill="firebrick4", alpha=0.2)+
  annotate(geom="text", label="Normalgewicht", color="darkgreen", x=65, y=22)+
  annotate(geom="text", label="Übergewicht", color="darkorange3", x=45, y=27.5)+
  annotate(geom="text", label="Adipositas (I)", color="firebrick2", x=45, y=32.5)+
  annotate(geom="text", label="Adipositas (II)", color="firebrick4", x=45, y=37.5)

P3
`geom_smooth()` using formula = 'y ~ x'

4 Modellprüfung

Auch hier prüfen wir die Annahmen bzw. Modellvoraussetzungen in grafischer Form:

par(mfrow=c(1,2))
plot(LM.BMI2, c(1:2))

Die Verteilung der Residuen erscheint zunächst etwas problematischer als im zwei-Gruppen Vergleich. Im Q-Q-Plot erkennen wir eine weitgehende Normalverteilung. Die Residuale aus dem ersten Plot weichen lediglich an den Rändern des Vorhersagebereichs deutlicher ab. Hier befinden sich allerdings auch eher wenig Datenpunkte – wichtiger ist der Bereich, in dem sich der Großteil der Daten befindet. Für den Großteil der Residuen wirkt die Verteilung dennoch trendfrei und zufällig, so dass wir auch in diesem Fall die wichtigsten Modellvoraussetzungen bzw. Annahmen als erfüllt ansehen können.

4.1 Ergebniszusammenfassung

Wir stellen einen statistisch signifikanten Zusammenhang zwischen dem Taillenumfang und dem BMI fest. Pro Prozent mehr Taillenumfang (relativ zur Körpergröße) nimmt der BMI im Schnitt um ca. 2.27 \(kg/m^2\) zu.

5 Fazit

Wir erkennen, dass t-Test und ALM, sowie einfache lineare Regression und ALM letztendlich jeweils austauschbare Verfahren sind und die Anwendung nur davon abhängt, welche Art von Prädiktorvariable wir einsetzen (Kategorisch/Gruppierung vs. metrisch).

Wir können sowohl den t-Test als auch die einfache lineare Regression durch das ALM ersetzen. Im letzteren Fall ist die Austauschbarkeit offensichtlich. Im Fall des t-Tests erschließt sich die Übereinstimmung mit dem ALM bei Beachtung der Codierung der Gruppenvariable mit den Werten 0 und 1, sowie den so ermittelten Koeffizienten \(b_0\) und \(b_1\). Wir werden in künftigen Einheiten sehen, dass das ALM noch weitere Verfahren in sich vereint, die wir bisher immer getrennt voneinander gehandhabt haben.


Ende | Stand: 2024-04-11