PDF-Version

library(readr)
library(pander)
library(broom)
library(dplyr)
library(ggplot2)

Aufgabe 1

Die Datei school.csv beinhaltet Daten zur Zulassung an amerikanischen Hochschulen. admit bezeichnet dabei, ob der jeweilige Aspirant an der Hochschule aufgenommen wurde, und gpa bezeichnet den Notendurchschnitt des Highschool-Abschlusses. Von Interesse ist, ob der Notendurchschnitt einen Einfluss auf die Aufnahmewahrscheinlichkeit an der Universität hat.

school <- read_csv("data/school.csv", col_types = cols()) %>%
  rename_all(tolower)

a)

Passen Sie ein lineares Regressionsmodell an die Daten an. Zeichnen Sie einen Scatterplot der Daten und die geschätzte Regressionsgerade.

mod1a <- lm(admit ~ gpa, data = school)
pander(mod1a)
Fitting linear model: admit ~ gpa
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.4224 0.2061 -2.05 0.04104
gpa 0.2183 0.06041 3.613 0.0003412
ggplot(school, aes(x = gpa, y = admit)) +
  geom_point() +
  geom_smooth(method = lm, color = "red", alpha = .1) +
  scale_x_continuous(breaks = seq(0, 4, .5), minor_breaks = seq(0, 4, .25)) +
  scale_y_continuous(breaks = 0:1, labels = c("No", "Yes")) +
  labs(
    title = "University admission by GPA",
    subtitle = "A terrible model.",
    x = "GPA", y = "Admitted"
  )

b)

Zeigen Sie anhand geeigneter Graphen der studentisierten Residuen, dass die Annahme homoskedastischer und normalverteilter Fehler verletzt ist.

mod1a %>%
  augment() %>%
  ggplot(aes(x = .fitted, y = .std.resid)) +
  geom_point() +
  labs(
    title = "University admission by GPA",
    subtitle = "Residuals + fitted values",
    x = "Fitted value", y = "Studentized Residual"
  ) -> plot_resid

mod1a %>%
  augment() %>%
  ggplot(aes(sample = .std.resid)) +
  geom_qq() +
  geom_qq_line() +
  coord_flip() +
  scale_x_continuous(breaks = seq(-4, 4, 1), minor_breaks = NULL) +
  scale_y_continuous(breaks = seq(-4, 4, 1), minor_breaks = NULL) +
  labs(
    title = "",
    subtitle = "QQ-plot of residuals",
    y = "Sample quantile", x = "Normal quantile"
  ) -> plot_qq

cowplot::plot_grid(plot_resid, plot_qq)

c)

Die Verletzung der Modellannahmen des linearen Modells kann leicht begründet werden. Nehmen Sie an, \(Y\) sei eine binäre Zufallsvariable mit den Ausprägungen \(0\) oder \(1\) und es gelte ein lineares Regressionsmodell, d. h. der Zusammenhang

\[Y = \beta_0 + \beta_1 x + \varepsilon\]

mit \(\mathbb{E}(\varepsilon) = 0\). Zeigen Sie, dass dann im Allgemeinen \(\mathrm{Var}(\varepsilon)\) nicht konstant ist, also hier von \(x\) abhängt.

Für eine normalverteilte Variable \(Y \sim \mathcal{N}(\mu, \sigma^2)\) gilt \(\mathbb{E}(Y) = \mu\) und \(\mathbb{V}(Y) = \sigma^2\), aber da in diesem Fall \(Y \sim \mathrm{Bernoulli}(p)\), gilt \(\mathbb{E}(Y) = p\) und \(\mathbb{V}(Y) = p (1 - p)\). Damit ist hier die Varianz abhängig vom Erwartungswert (dessen Schätzung wiederum abhängig von den Daten ist), im Gegensatz zu einer normalverteilten Variable, in der Erwartungswert und Varianz unabhängig sind.

\[\begin{aligned} \mathbb{E}(Y) &= \beta_0 + \beta_1 x = p && 0 \leq p \leq 1 \\\\ V(Y) &= V(\beta_0 + \beta_1 x + \varepsilon) \\ p (1-p) &= V(\beta_0 + \beta_1 x + \varepsilon) && \bigg|\ \text{da } Y \sim \mathrm{Bernoulli}(p) \\ p (1-p) &= V(\varepsilon) \end{aligned}\]

Aufgabe 2

Zur Bestimmung der Toxizität des Digitalispräparats Ouabain bei Fröschen wurden für verschiedene Dosierungen \(z\) (in µg pro 100 g Froschgewicht) jeweils \(n = 30\) Frösche mit der Dosis \(z\) geimpft, und dann die Anzahl \(r\) der am Gift gestorbenen Tiere gezählt:

quabain <- tibble(
  z = c(23, 25, 27, 29, 31, 33, 35, 37),
  r = c(0, 1, 2, 7, 18, 25, 28, 29)
)
Ouabain: Dosis und Todesfälle
z 23 25 27 29 31 33 35 37
r 0 1 2 7 18 25 28 29

Die beobachtete Toxizitätsrate ist dann die relative Häufigkeit \(y = r / n\). Die Abhängigkeit der erwarteten Toxizitätsrate \(\mu(x) = \mathbb{E}(Y\ |\ x)\) von der log-Dosis \(x = \log_{10}(z)\) soll durch ein logistisches Modell beschrieben werden:

\[ \operatorname{logit} \mu(x) = \beta_1 + \beta_2 x = \eta(x) \]

Schätzen Sie den Parameter \(\beta\) nach den folgenden beiden Methoden:

a)

Bestimmen Sie die \((x, v)\)-Regressionsgerade für die log-Dosen \(x_i\) und die ungewichteten transformierten Beobachtungen \(v_i = \operatorname{logit}(y_i)\), aber nur für \(i \geq 2\), also ohne \(v_1 = \operatorname{logit}(0) = -\infty\).

quabain_a2 <- quabain %>%
  mutate(
    x = log10(z),         # log-Dosis
    y = r / 30,           # Anteilswert
    v = log(y / (1 - y)), # Logit-Anteilswert
    w = 30 * y * (1 - y)  # Gewicht
  ) %>%
  filter(is.finite(v)) # Entfernt infinite Werte in v, i.e. für y = 0

mod2a <- lm(v ~ x, data = quabain_a2)
pander(mod2a, caption = "Quabain, ungewichtetes Modell")
Quabain, ungewichtetes Modell
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -62.85 2.77 -22.69 3.09e-06
x 42.32 1.86 22.75 3.05e-06
ggplot(quabain_a2, aes(x = x, y = v)) +
  geom_point() +
  geom_smooth(method = lm, color = "red", alpha = .1) +
  labs(
    title = "Quabain", subtitle = "Logit/Log Modell, ungewichtet",
    x = "Log-Dosis (x)", y = "Logit-proportion (v)"
  )

b)

Bestimmen Sie die gewichtete \((x,v)\)-Regressionsgerade für die gewichteten Punkte \((x_i, v_i, w_i)\) mit Gewichten \(w_i = n_i y_i (1 - y_i)\), aber nur für \(i \geq 2\), weil \(w_1 = 0\) (Das Gewicht \(w_i\) ist eine Approximation für die reziproke Varianz der transformieren Zufallsvariablen \(V_i = \operatorname{logit}(Y_i)\)).

mod2b <- lm(v ~ x, weights = w, data = quabain_a2)
pander(mod2b, caption = "Quabain, gewichtetes Modell")
Quabain, gewichtetes Modell
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -66.12 3.304 -20.01 5.763e-06
x 44.53 2.22 20.06 5.692e-06
ggplot(quabain_a2, aes(x = x, y = v)) +
  geom_point() +
  geom_smooth(method = lm, color = "red", alpha = .1, aes(weight = w)) +
  labs(
    title = "Quabain", subtitle = "Logit/Log Modell, gewichtet",
    x = "Log-Dosis (x)", y = "Logit-proportion (v)"
  )

Vergleichen Sie die beiden Schätzungen auch graphisch. Plotten Sie dazu die geschätzen Gerade \(\hat{\eta}(x)\) in die transformierten Punkte \((x, v)\), sowie die geschätzten Funktionen \(\hat{\mu}(x)\) in die Punkte \((x, y)\).

Notiz:

\[\begin{aligned} \operatorname{logit}(x) &= \log\left(\frac{x}{1-x}\right)\\ \operatorname{logit}^{-1}(x) = \operatorname{logistic}(x) &= \frac{1}{1 + \exp(-x)} \end{aligned}\]

both_models <- bind_rows(
  augment(mod2a) %>% mutate(Modell = "Ungewichtet"),
  augment(mod2b) %>% mutate(Modell = "Gewichtet")
)

both_models %>%
  ggplot(aes(x = x, y = v)) +
  geom_point() +
  geom_line(aes(y = .fitted, color = Modell)) +
  scale_color_brewer(palette = "Set1", guide = FALSE) +
  labs(
    title = "Modellvergleich",
    subtitle = "Auf Skala des linearen Prädiktors",
    x = "Log-Dosis (x)", y = "Logit-proportion (v)"
  ) -> plot_linear

both_models %>%
  mutate(.fitted_transf = 1 / (1 + exp(-.fitted))) %>%
  left_join(quabain_a2, by = c("x", "v")) %>%
  ggplot(aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = .fitted_transf, color = Modell)) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "",
    subtitle = "Auf Skala der Response",
    x = "Log-Dosis (x)", y = "Proportion (y)"
  ) -> plot_trans

cowplot::plot_grid(plot_linear, plot_trans)

Aufgabe 3

Zeigen Sie, dass im allgemeinen linearen Modell für den ungewichteten Kleinste-Quadrate-Schätzer \(\hat{\beta}\) von \(\beta\) die folgenden Aussagen gelten:

\[\begin{aligned} \mathbb{E}(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta} && \mathrm{Cov}(\hat{\boldsymbol{\beta}}) = \sigma^2 (\boldsymbol{X}^T\boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{W} \boldsymbol{X} (\boldsymbol{X}^T\boldsymbol{X})^{-1} \end{aligned}\]

\[\begin{aligned} \mathbb{E}(\widehat{\boldsymbol{\beta}}) &= \mathbb{E}\left( \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \boldsymbol{Y} \right) \\ &= \mathbb{E}\left( \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \left(\boldsymbol{X} \boldsymbol{\beta}+ \epsilon \right) \right) \\ &= \mathbb{E}\left( \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{\beta} + \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \boldsymbol{\varepsilon} \right) \\ &= \underbrace{\left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \boldsymbol{X}}_{=\ \boldsymbol{I}} \boldsymbol{\beta} + \underbrace{\left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \mathbb{E}(\boldsymbol{\varepsilon})}_{=\ 0} & \bigg\vert\ \mathbb{E}(\boldsymbol{\varepsilon}) = 0\\ &= \boldsymbol{\beta} \\&& \Box \end{aligned}\]

\[\begin{aligned} \mathrm{Cov}(\widehat{\boldsymbol{\beta}}) &= \mathrm{Cov}\left( \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \boldsymbol{Y} \right) \\ &= \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \mathbb{V}(\boldsymbol{Y}) \left(\left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T\right)^T \\ &= \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \mathbb{V}(\boldsymbol{Y}) \boldsymbol{X} {\left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1}}^T \\ &= \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \mathbb{V}(\boldsymbol{X}\boldsymbol{\beta} + \varepsilon) \boldsymbol{X} {\left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1}}^T \\ &= \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \mathbb{V}(\varepsilon) \boldsymbol{X} {\left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1}}^T \\ &= \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \sigma^2\boldsymbol{W} \boldsymbol{X} \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \\ &= \sigma^2 \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \boldsymbol{X}^T \boldsymbol{W} \boldsymbol{X} \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1} \\ && \Box \end{aligned}\]