8 Lineare Modelle

Warning

Hinweis: Achtung, die Ergebnisse können selbst mit set.seed() etwas variieren!

8.1 Aufgabe 1

Das folgende Modell soll die Größe der Menschen einer Population anhand des Körpergewichts vorher sagen. Welches der folgenden Modelle ist am besten für den genannten Zweck definiert?

a)
\[ \begin{aligned} height_{i}\sim Normal(\mu_{i},\sigma) \\ \mu_{i} = \alpha + \beta \cdot weight_{i} \\ \alpha \sim Normal(178, 20) \\ \beta \sim Normal(0,10) \\ \sigma \sim Exp(0.1) \end{aligned} \]

b)
\[ \begin{aligned} height_{i}\sim Normal(\mu,\sigma)\\ \mu_{i} = \alpha + \beta \cdot weight\\ \alpha \sim Normal(178, 20)\\ \beta \sim Normal(0,10)\\ \sigma \sim Exp(0.1) \end{aligned} \]

c)
\[ \begin{aligned} height_{i}\sim Normal(\mu_{i},\sigma)\\ \mu_{i} = \alpha + \beta \cdot weight_{i}\\ \alpha \sim Normal(178, 20)\\ \beta \sim Normal(5,3)\\ \sigma \sim Exp(0.1) \end{aligned} \]

d)
\[ \begin{aligned} height_{i}\sim Normal(\mu_{i},\sigma)\\ \mu_{i} = \alpha + \beta \cdot weight_{i}\\ \alpha \sim Normal(178, 20)\\ \beta \sim Normal(40,100)\\ \sigma \sim Exp(0.1) \end{aligned} \]

Antwort A:
Die Modelldefinition der Option A ist zwar grundsätzlich nicht falsch, geht aber davon aus, dass der Anstieg unseres linearen Modells \(\beta \sim N(0,10)\) positiv sowie negativ sein könnte. Genauer gesagt, gibt die angegebene Normalverteilung an, dass die Wahrscheinlichkeit gleich groß für positive und negative Werte ist (\(\mu = 0,\pm1\ SD = -10 bis +10\)). Praktisch interpretiert heißt das: Das Modell nimmt an, dass Menschen, die immer schwerer werden, mit gleich großer Wahrscheinlichkeit größer oder klein er werden könnten. Eine bessere Modelldefinition wäre allerdings, dass Menschen, die immer schwerer werden auch mit größerer Wahrscheinlichkeit größer werden.

Antwort B:
Modell B enthält einen Formfehler. Sowohl im Likelihood als auch im linearen Modell fehlen die Beobachtungen \(i\). Entsprechend wird laut Modelldefinition keine konkrete Beobachtung vorhergesagt, was ja unter anderem der Sinn eines Modells ist.

Antwort C:
Antwort C ist richtig.
Während Antwort A und C zwar beide nicht vollkommen falsch sind, kann man im Anstieg des Modells \(\beta\) erkennen, dass dieser apriori überwiegend positiv definiert ist. Der Priori-Anstieg geht also davon aus, dass die Größe einer Person zunimmt, je schwerer sie ist. Als Modell, das die Größe einer Person anhand des Gewichts vorhersagen soll, ist dieses von den aufgelisteten am besten geeignet.

Antwort D:
In Modell D ist der Anstieg \(\beta\) unrealistisch groß definiert. Das Modell nimmt apriori an, dass es möglich ist, dass ein Mensch mit einem Gewicht von \(x=0\) sich von einem Menschen mit \(x=1\) um 140 cm unterscheidet. Praktisch gesprochen:
Wenn wir die Körpergröße zweier Menschen vergleichen, die sich in ihrem Körpergewicht nur um einen Kilogramm unterscheiden (Mensch 1 = 60 kg, Mensch 2 = 61 kg), dann wäre es laut dem Modell möglich, dass Mensch 2 auf Grundlage seines höheren Gewichts 140 cm größer sein könnte als Mensch 1. Da das in der Realität aber (fast) nie der Fall sein wird, handelt es sich um keine plausible Annahme.

8.2 Aufgabe 2

Laden Sie für die folgenden Aufgaben den Datensatz “kidiq” aus dem Paket rstanarm.

data(kidiq)

Erstellen Sie ein Modell, dass den IQ von Kindern anhand des IQs der Mutter vorhersagt. Was sagt der Intercept des Modells aus?

set.seed(42)
m1a <- stan_glm(kid_score ~ mom_iq, data = kidiq, refresh = 0) #Modell formulieren
parameters(m1a)

Das Modell sagt den IQ der Kinder anhand des IQs der Mutter voraus. Der Intercept des Modells stellt wie immer den erwarteten (mittleren) y-Wert für eine Beobachtung mit x = 0 dar. Das heißt wir sehen den erwarteten mittleren IQ eines Kindes, von einer Mutter, die einen IQ von 0 hat. Da das in der Realität nicht möglich ist, ist das Modell nur begrenzt gut interpretierbar.

8.2.2 Aufgabe 2 b)

A2 b)
L2 b)

Zentrieren Sie die Prädiktorvariable (x-Variable/UV) und stellen Sie ein Modell auf, das den IQ von Kindern mit zentrierten Werten vorhersagt. Wie ist der Intercept nun zu interpretieren?

kidiq <- 
  kidiq%>%
  mutate(mom_iq_c = mom_iq - mean(mom_iq)) #Ein Wert in der Spalte mom_iq minus den Mittelwert der Spalte mom_iq

set.seed(42)
m1b <- stan_glm(kid_score ~ mom_iq_c, data = kidiq, refresh = 0) #Modell mit zentriertem Prädiktor
parameters(m1b)

Jetzt, da die Prädiktorvariable zentriert ist, können wir den Intercept besser verstehen. Nun sagt dieser aus, dass das Modell im Median einen mittleren IQ von 86,8 bei einem Kind erwartet, wenn die Mutter einen durchschnittlichen IQ hat.

8.2.3 Aufgabe 2 c)

A2 c)
L2 c)

Innerhalb welcher beiden Werte liegt das 90 % PI des mittleren IQs der Kinder?

m1btab <- m1b%>% #Modell m1b als Tabelle ausgeben
  as_tibble()

m1btab%>%
  select(`(Intercept)`)%>%
  eti(ci = .9)

8.2.4 Aufgabe 2 d)

A2 d)
L2 d)

Welcher Wert bildet die Unsicherheit hinsichtlich des mittleren IQs der Kinder am Intercept ab?

m1b

stan_glm
 family:       gaussian [identity]
 formula:      kid_score ~ mom_iq_c
 observations: 434
 predictors:   2
------
            Median MAD_SD
(Intercept) 86.8    0.8  
mom_iq_c     0.6    0.1  

Auxiliary parameter(s):
      Median MAD_SD
sigma 18.3    0.6  

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Wenn man sich den Output des Modells anguckt, sieht man, dass der mittlere IQ der Kinder am Intercept einen Median von 86.8 hat. Allerdings hat die Verteilung auch eine mediane Standardabweichung (MAD_SD) von 0.9. Der Wert der die Unsicherheit zum mittleren IQ der Kinder am Intercept darstellt, ist also 0.8.

8.2.5 Aufgabe 2 e)

A2 e)
L2 e)

Wie groß ist die Wahrscheinlichkeit, dass der Regressionskoeffizient größer ist als 0,65?

m1btab%>%
  summarise(mean(mom_iq_c > .65))

Die Wahrscheinlichkeit, dass der Anstieg größer ist als 0.65 liegt bei ca. 25 Prozent.

8.2.6 Aufgabe 2 f)

A2 f)
L2 f)

Wie groß ist die Wahrscheinlichkeit, dass der mittlere IQ eines Kindes, dessen Mutter fünf Einheiten über dem Durchschnitt liegt, größer ist als 90?

m1btab_f <-
  m1btab%>%
  mutate(x5 = `(Intercept)` + mom_iq_c*5) #Spalte mutieren, in der man zum Intercept (x = 0), fünf Einheiten von x hinzufügt

m1btab_f%>%
  count(x5 > 90)%>%
  mutate(Anteil = n / sum(n))

Die Wahrscheinlichkeit, dass der mittlere IQ eines Kindes, das von einer Mutter mit x = 5 stammt, größer ist als 90, liegt bei ca. 42 Prozent.

8.3 Aufgabe 3

Laden Sie für die folgenden Aufgaben den Datensatz “mtcars” aus dem Paket tidyverse.

data(mtcars)

8.3.1 Aufgabe 3 a)

A3 a)
L3 a)

Erstellen Sie ein Modell, dessen AV der Spriverbrauch und dessen UV das zentrierte Gewicht eines Autos ist.

mtcars <-
  mtcars%>%
  mutate(weight_c = center(wt)) #Gewicht zentrieren

set.seed(42)
m2a <-
  stan_glm(mpg ~ weight_c, data = mtcars, refresh = 0)
m2a

stan_glm
 family:       gaussian [identity]
 formula:      mpg ~ weight_c
 observations: 32
 predictors:   2
------
            Median MAD_SD
(Intercept) 20.1    0.5  
weight_c    -5.3    0.5  

Auxiliary parameter(s):
      Median MAD_SD
sigma 3.1    0.4   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

8.3.2 Aufgabe 3 b)

A3 b)
L3 b)

Wie groß ist die Wahrscheinlichkeit, dass der mittlere Spritverbauch kleiner ist als 20 mpg?

m2atab <- 
  m2a%>% #Modell als Tabelle
  as_tibble()

m2atab%>%
  summarise(mean(`(Intercept)` < 20))

Die Wahrscheinlichkeit, dass der mittlere Spritverbrauch kleiner ist als 20 mpg, liegt bei 43 Prozent.

8.3.3 Aufgabe 3 c)

A3 c)
L3 b)

Wie groß ist die Wahrscheinlichkeit, dass eine mittlere Streuung von 3,5 nicht überschritten wird?

m2atab%>%
  summarise(mean(sigma < 3.5))

Die Wahrscheinlichkeit, dass eine mittlere Streuung von 3,5 nicht überschritten wird, liegt bei ca. 82 Prozent.

8.4 Codesammlung

Tipps für eine Codesammlung

Parameter eines Modells anzeigen:

parameters(m1)

UV zentrieren:

d <- 
  d%>%
  mutate(uv_c = uv - mean(uv))

X % PI angeben

# Wichtig: Erst in Tabelle umwandeln!
m1tab <- m1%>%
  as_tibble()

m1tab%>%
  select(`(Intercept)`)%>%
  eti(ci = .X)

Wahrscheinlichkeit, dass Regressionskoeffizient größer als X ist:

m1tab%>%
  summarise(mean(uv_c > X))

Wahrscheinlichkeit, dass der mittlere Wert der AV, wenn die UV X Einheiten über dem Durchschnitt liegt, größer als Y ist

m1tab_f <-
  m1btab%>%
  mutate(xX = `(Intercept)` + uv_c*5) 

m1tab_f%>%
  count(xX > Y)%>%
  mutate(Anteil = n / sum(n))

Wahrscheinlichkeit, dass die mittlere Streuung von X nicht überschritten wird

m1tab%>%
  summarise(mean(sigma < X))