PDF-Version

Aufgabe 1

Ziel einer retrospektiven epidemiologischen Studie war es, den prophylaktischen Effekt von Weinkonsum auf eine bestimmte Herzerkrankung zu untersuchen (0 = “keine Herzerkrankung”, 1 = “Herzerkrankung”). Des Weiteren wurde beobachtet, ob die Patienten in den letzten zehn Jahren geraucht haben oder nicht. Eine logistische Regressionsanalyse mit zwei binären Einflussgrößen “Weinkonsum” und “Rauchen” (jeweils 0 = Nein, 1 = Ja) und dem Interaktionsterm Weinkonsum×Rauchen ergab folgenden Output

Variable Schätzer Odds Ratio
Intercept -3.3
Weinkonsum -0.3 0.74
Rauchen 0.6 1.82
Weinkonsum × Rauchen 1.2 3.32

a)

Interpretieren Sie die Schätzungen der Odds Ratios.

Das Modell hat erstmal die Form

\[\begin{aligned} \log\left(\frac{\mathrm{P}(\mathrm{Erkrankt})}{1 - \mathrm{P}(\mathrm{Erkrankt})}\right) &= -3.3 - 0.3 \cdot \mathrm{Weinkonsum} + 0.6 \cdot \mathrm{Rauchen} + 1.2 \cdot \mathrm{Weinkonsum} \cdot \mathrm{Rauchen} \\ \frac{\mathrm{P}(\mathrm{Erkrankt})}{1 - \mathrm{P}(\mathrm{Erkrankt})} &= \exp(-3.3) \exp(-3.3 \cdot \mathrm{Wein}) \exp(0.6 \cdot \mathrm{Rauchen}) \exp(1.2 \cdot \mathrm{Wein} \cdot \mathrm{Rauchen}) \end{aligned}\]

  • Intercept: Wird nicht mitinterpretiert
  • Weinkonsum: Weinkonsumierende Raucher haben ein um den Faktor \(0.74 \cdot 3.32 \approx 2.46\) erhöhtes Risiko verglichen mit Nichttrinkenden Nichtrauchern, und weintrinkende Nichtraucher ein um den Faktor \(0.74\) verringertes Risiko.
  • Rauchen: Rauchende Nichttrinker haben ein um den Faktor 1.82 erhöhtes Risiko, Rauchende Weintrinker wiederum ein um \(1.82 \cdot 3.32 \approx 6.05\) höheres Risiko als nichtrauchende Trinker.
  • Weinkonsum×Rauchen: (Multiplikative) Veränderung des OR für weintrinkende Raucher verglichen mit weintrinkenden Nichtrauchen, nichttrinkenden Rauchern oder nichttrinkenden Nichtrauchern.

(Anm. v. Lukas: Okay, jetzt hab ich einen Knoten im Kopf.)

b)

Wie lautet das geschätzte Odds Ratio von weintrinkenden Rauchern im Vergleich zu nichttrinkenden Nichtrauchern, und wie lautet das geschätzte Odds Ratio von weintrinkenden Rauchern zu weintrinkenden Nichtrauchern?

Für den Vergleich von weintrinkenden Rauchern und nichttrinkenden Nichtrauchern reicht es nicht, einen einzelne Koeffizienten zu betrachten, und da wir zusätzlich den Interaktionsterm im Modell haben, müssen wir die drei Terme zusammenfassen:

\[\begin{aligned} \log\left(\frac{\mathrm{P}(\mathrm{Erkrankt})}{1 - \mathrm{P}(\mathrm{Erkrankt})}\right) &= -3.3 - 0.3 + 0.6 + 1.2 \\ \frac{\mathrm{P}(\mathrm{Erkrankt})}{1 - \mathrm{P}(\mathrm{Erkrankt})} &= \exp(-3.3) \exp(-0.3) \exp(0.6) \exp(1.2)\\ \frac{\mathrm{P}(\mathrm{Erkrankt})}{1 - \mathrm{P}(\mathrm{Erkrankt})} &= \exp(-3.3) \cdot 4.4817 \end{aligned}\]

Damit entspricht das Odds Ratio für die weintrinkenden Nichtraucher \(\approx 4.5\), also einmal etwa 4.5fach höherem Risiko als die nichttrinkenden Nichtraucher.

Für den Vergleich der weintrinkenden Rauchern zu weintrinkenden Nichtrauchern gehen wir ähnlich vor, allerdings der einfacheren Betrachtung halber in zwei Schritten: Wir berechnen die (log)-Odds für beide uns interessierenden Gruppen (Weintrinken = 1 & Rauchen = 1 vs Weintrinken = 1 & Rauchen = 0), teilen die Odds und erhalten damit as OR:

\[\begin{aligned} \log(\mathrm{Odds})_{\text{Wein+Rauchen}} &= -3.3 - 0.3 + 0.6 + 1.2 \\ \log(\mathrm{Odds})_{\text{Wein+Nichtrauchen}} &= -3.3 - 0.3 \\ \mathrm{Odds Ratio} &= \frac{\exp(-3.3) \exp(-0.3) \exp(0.6) \exp(1.2)}{\exp(-3.3) \exp(-0.3)} \\ &= \exp(0.6)\exp(1.2) \\ &\approx \underline{\underline{6.05}} \end{aligned}\]

(Anm. v. Lukas: Das war nicht so mehrschrittig nötig, aber der Knoten in meinem Kopf tut jetzt weniger weh.)

c)

Geben Sie für alle möglichen Kovariablenkombinationen die geschätzte Wahrscheinlichkeit an, an dem Herzleiden zu erkranken. Wie ist die Aussage zu beurteilen, die Gefahr einer Erkrankung an dem entsprechenden Herzleiden lasse sich durch regelmäßigen Weinkonsum reduzieren?

Erstmal die Tabelle für die Log-Odds (\(\eta\)), im zweiten Schritt dann auf Skala der Response.

Kein Weintrinker Weintrinker
Nichtraucher -3.3 -3.3 - 0.3 = -3.6
Raucher -3.3 + 0.6 = -2.7 -3.3 - 0.3 + 0.6 + 1.2 = -1.8

Nun durch Anwendung der logistischen Funktion auf die jeweiligen \(\eta_{ij}\):

Kein Weintrinker Weintrinker
Nichtraucher 0.036 0.027
Raucher 0.067 0.142

Das Erkrankungsrisiko lässt sich also wohl nicht durch Weinkonsum an sich reduzieren, viel eher sehen wir, dass die Kombination aus Nichtrauchen und Weinkonsum die geringste, während Rauchen mit Weinkonsum die größte Erkrankungswahrscheinlichkeit aufweist. Weinkonsum wäre in dem Sinne also nur risikomindernd, wenn man dabei auch nicht raucht oder Rauchentwöhnung gleich mitempfiehlt.

d)

Wie wären wohl die Parameterschätzungen ausgefallen, hätte man die Interaktion nicht berücksichtigt? Begründen Sie Ihre Vermutung.

Wir hätten vermutlich einen größeren Koeffizienten für Weinkonsum erhalten, da wir das erhöhte Risiko durch Rauchen sozusagen mit aufgefangen hätten. Die Koeffizienten müssten sich ja aus dem obigen Modell mit Interaktionsterm ableiten lassen.

Aufgabe 2

a)

Wir betrachten den Fall der binären Regression mit gruppierten Daten. Überzeugen Sie sich von der in der Vorlesung getätigten Aussage, dass für den Fall von stochastisch unabhängigen \(Y_l \sim B(n_l, \pi_l)\) für \(l = 1, \ldots, m\) gilt, dass die log-Likelihood von \((Y_1, \ldots, Y_m)\) bis auf konstante Terme, die von \(\pi_l\) unabhängig sind, die folgende Form hat:

\[ \sum_{l=1}^m \left\{Y_l \log(\pi_l) - Y_l \log(1-\pi_l) + n_l \log(1 - \pi_l)\right\} \]

\[\begin{aligned} f(Y_i\ |\ \boldsymbol{x}_i \beta_i) &= \binom{n}{Y_i} \pi_i^{Y_i} (1 - \pi_i)^{n - Y_i}\\[1em] L(\beta) &= \prod^n_{i=1} f(Y_i\ |\ \boldsymbol{x}_i \beta_i) \\ \Longrightarrow l(\beta) &= \sum^n_{i=1} \log f(Y_i\ |\ \boldsymbol{x}_i \beta_i) \\ &= \sum^n_{i=1} \log\left(\binom{n}{Y_i} \pi_i^{Y_i} (1 - \pi_i)^{n - Y_i}\right) \\ &= \sum^n_{i=1} \log\left(\binom{n}{Y_i}\right) + Y_i \cdot \log(\pi_i) + (n - Y_i) \cdot \log(1 - \pi_i) \\ &= \sum^n_{i=1} \log\left(\binom{n}{Y_i}\right) + Y_i \cdot \log(\pi_i) - Y_i \cdot \log(1 - \pi_i) + n \cdot \log(1 - \pi_i) \\\\&&\Box \end{aligned}\]

b)

Zeigen Sie, dass für die Linkfunktion \(g(\pi) = \log(\frac{\pi}{1 - \pi})\) des Logit-Modells

\[ g(1-\pi) = -g(\pi) \]

gilt.

Durch Einsetzen in die logistische Funktion erhalten wir jeweils den selben Term, im zweiten Teil durch Anwendung der Quotientenregel des Logarithmus:

\[\begin{aligned} g(1 - \pi) &= \log\left(\frac{1 -\pi}{1 - (1 - \pi)}\right) &= \log\left(\frac{1 -\pi}{\pi}\right) \\ -g(\pi) &= - \log\left(\frac{\pi}{1 - \pi}\right) &= \log\left(\frac{1 - \pi}{\pi}\right) \end{aligned}\]

c)

Zeigen Sie, dass für die Linkfunktion \(g(\pi) = -\log(-\log(\pi))\) des Log-Log-Modells

\[ g(\pi) \neq -g(1 - \pi) \]

gilt.

\[\begin{aligned} \text{Sei } \pi = 0.2 \\ -g(0.2) &= \log(-\log(0.2)) \approx 0.476 \\ g(1 - 0.2) &= -\log(-\log(0.8)) \approx 1.5 \end{aligned}\]

Aufgabe 3

Die Datei school2.sas beinhaltet Daten zur Zulassung an amerikanischen Hochschulen. admit bezeichnet dabei, ob der jeweilige Aspirant an der Hochschule aufgenommen wurde, gpa bezeichnet den Notendurchschnitt des Highschool-Abschlusses, gre den Wert eines Aufnahmetests (Skala: 200 bis 800) und rank das Renommee der Highschool des Aspiranten (1=exzellentes Renommee, 2=gutes Renommee, 3=mittleres Renommee und 4=schlechtes Renommee). Von Interesse ist die Wahrscheinlichkeit einer Hochschulzulassung.

a)

Passen Sie zunächst ein logistisches Modell an die Daten mit den Einflussvariablen gre, gpa und rank an die Daten an. Interpretieren Sie die Effektschätzer.

library(readr)
library(broom)

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

mod2a <- glm(admit ~ gre + gpa + rank, family = binomial(), data = school2)

mod2a %>%
  tidy() %>%
  mutate(OR = exp(estimate)) %>%
  select(Variable = term, Schätzer = estimate, OR) %>%
  kable2() %>%
  kable_styling2()
Variable Schätzer OR
(Intercept) -3.450 0.032
gre 0.002 1.002
gpa 0.777 2.175
rank -0.560 0.571
  • gre: Da das OR sehr nah an 1 ist scheint das Ergebnis des Aufnahmetests kaum zu einem Unterschied der Aufnahmewahrscheinlichkeitzu führen, wenn gpa und rank gleich snd.
  • gpa: Aspiranten mit einem um 1 höheren GPA als andere Aspiranten haben eine 2.1-fach höhere Aufnahmewahrscheinlichkeit (wieder bei gleichem Aufnahmetestergebnis und Renommee der Hochschule).
  • rank: Highschools mit einem um 1 höheren Rang / niedrigerem Renommee-Grad zeigen eine etwa halb so große Aufnahmewahrscheinlichkeit wie andere Hochschulen, analog wieder bei gleichem GPA und Aufnahmetestergebnis.

b)

Passen Sie ein logistisches Modell an die Daten an, aber die Variable rank soll über die Dummy-codierten Variablen rank1, rank2 und rank3 ins Modell eingehen (d.h. rank1 = 1 \(\leftrightarrow\) rank = 1, rank2 = 1 \(\leftrightarrow\) rank = 2, rank3 = 1 \(\leftrightarrow\) rank = 3).

Interpretieren Sie die Effektschätzer zu den Dummy-Variablen. Finden Sie das Modell aus a) oder b) geeigneter?

school2 <- school2 %>%
  mutate(rank_cat = relevel(factor(rank), ref = 4))

mod2b <- glm(admit ~ gre + gpa + rank_cat, family = binomial(), data = school2)

mod2b %>%
  tidy() %>%
  mutate(OR = exp(estimate)) %>%
  select(Variable = term, Schätzer = estimate, OR) %>%
  kable2() %>%
  kable_styling2()
Variable Schätzer OR
(Intercept) -5.541 0.004
gre 0.002 1.002
gpa 0.804 2.235
rank_cat1 1.551 4.718
rank_cat2 0.876 2.401
rank_cat3 0.211 1.235

Da es sich um 4 Kategorien handelt, deren numerische Rang wohl eher nicht wie eine metrische Variable interpretierbar ist, wird so der Vergleich der verschiedenen Rangstufen einfacher. In diesem Modell werden die ORs für die Renommee-Ränge 1-3 jeweils in Relation zu Rang 4 (also dem schlechtesten Renommee) interpretiert. Daher sehen wir, dass Aspiranten von sehr angesehenen Highschools eine 4 mal so große Aufnahmewahrscheinlichkeit als Aspiranten von den am wenigsten angesehen Highschools haben.

c)

Führen Sie mit dem Modell aus b) die folgende lineare Hypothesentests mit einem Signifikanzniveau von 5% durch (TEST oder CONTRAST Statement in SAS):

  1. Keine der Variablen hat einen Einfluss auf die Wahrscheinlichkeit für eine Hochschulzulassung.

Mittels Likelihood-Ratio Test: “Volles” Modell gegen Null-Modell:

mod0 <- glm(admit ~ 1, family = binomial(), data = school2)

anova(mod0, mod2b, test = "LRT") %>% 
  pander::pander(caption = "LRT: Modell gegen Null-Modell")
LRT: Modell gegen Null-Modell
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
399 500 NA NA NA
394 458.5 5 41.46 7.578e-08

Das signifikante Ergebnis sagt uns, dass wir die Hypothese (Keine Variable hat einen Einfluss auf die Zulassung) verwerfen können.

  1. Ein exzellentes Renommee der Highschool hat einen gleich großen Einfluss wie ein gutes Renommee auf die Wahrscheinlichkeit für eine Hochschulzulassung.

Mittels Wald-Test:

# Kontrastmatrix: Koeffizienten 4 - 5
contrast_matrix <- matrix(c(0, 0, 0, 1, -1, 0), nrow = 1)
aod::wald.test(b = coef(mod2b), Sigma = vcov(mod2b), L = contrast_matrix)
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 4.6, df = 1, P(> X2) = 0.033

Die beiden Rangstufen unterscheiden sich also signifikant (\(p < 0.05\)) bzgl. der Hochschulzulassungswahrscheinlichkeit.

Aufgabe 4

Das einfachste logistische Modell modelliert den Einfluss einer binären Variable \(X\) auf die binäre Variable \(Y\). \(X\) kann dabei z. B. eine dichotome Expositionsvariable darstellen (0 = nicht exponiert, 1 = exponiert) und \(Y\) der Indikator für das Vorliegen einer Krankheit (0 = nicht erkrankt, 1 = erkrankt) sein. Das zugehörige logistische Modell ist

\[ \mathrm{Pr}(Y = 1\ |\ X = x) = \frac{\exp(\beta_0 + x \beta_1)}{1 + \exp(\beta_0 + x \beta_1)} \]

Modell (1) kann auch über eine 2×2-Tabelle dargestellt werden. Dabei sind in den einzelnen Zellen die Anteile an Erkrankten bzw. nicht Erkrankten dargestellt, die man auf Grund des Modells bei Exponierten und nicht Exponierten erwarten würde:

Y
1 0 Summe
\(X\) 1 \(\frac{\exp(\beta_0 + x \beta_1)}{1 + \exp(\beta_0 + x \beta_1)}\) \(\frac{1}{1 + \exp(\beta_0 + x \beta_1)}\) 1
0 \(\frac{\exp(\beta_0)}{1 + \exp(\beta_0)}\) \(\frac{1}{1 + \exp(\beta_0)}\) 1

In einer Studie kennt man \(\beta_0\) und \(\beta_1\) nicht und muss diese aus einer beobachteten Tafel mit Häufigkeiten Schätzen:

Y
1 0 Summe
\(X\) 1 \(a\) \(b\) \(a + b\)
0 \(c\) \(d\) \(c + d\)

1)

Zeigen Sie, dass in einer solchen Situation der Standard-Schätzer \(\log\left(\tfrac{ad}{bc}\right)\) für das Log-Odds-Ratio genau dem Maximum-Likelihood-Schätzer für \(\beta_1\) entspricht. Berechnen Sie auch den Maximum-Likelihood-Schätzer für \(\beta_0\).

\(\log\left(\frac{ad}{bc}\right)\) reduziert zu \(\beta_1\) für \(x = 1\):

\[\begin{aligned} \frac{a}{b} &= \frac{\frac{\exp(\beta_0 + x \beta_1)}{1 + \exp(\beta_0 + x \beta_1)}}{\frac{1}{1 + \exp(\beta_0 + x \beta_1)}} = \frac{\exp(\beta_0 + x \beta_1) \cdot (1 + \exp(\beta_0 + x \beta_1))}{1 + \exp(\beta_0 + x \beta_1)} = \exp(\beta_0 + x \beta_1) \\ \frac{c}{d} &= \frac{\frac{\exp(\beta_0)}{1 + \exp(\beta_0)}}{\frac{1}{1 + \exp(\beta_0)}} = \frac{\exp(\beta_0) \cdot (1 + \exp(\beta_0))}{1 + \exp(\beta_0)} = \exp(\beta_0) \end{aligned}\]

\[\begin{aligned} \log\left(\frac{ad}{bc}\right) &= \log\left(\frac{\frac{a}{b}}{\frac{c}{d}}\right) \\ &= \log\left(\frac{\exp(\beta_0 + x \beta_1)}{\exp(\beta_0)}\right) \\ &= \log\left(\exp(x \beta_1)\right) \\ &= x \beta_1 \end{aligned}\]

…allerdings hilft uns das nicht wirklich weiter. Die Bestimmung des MLE haben wir leider nicht hinbekommen, weil wir entweder am Nullsetzen der Score-Funktion gescheitert sind (Lukas) oder der Plan, den Ausdruck in die Form von \(\log\left(\frac{ad}{bc}\right)\) zu pressen (Roman), nicht aufging.

2)

Berechnen Sie die beobachtete Fisher-Information und daraus einen Varianzschätzer für den ML-Schätzer.

Hinweis: Die beobachtete Fisher-Matrix ist \(H(\beta) = - \frac{\partial^2}{\partial \beta \partial \beta'} l(\beta)\), wobei \(l(\beta)\) die Likelihood darstellt. Ein geeigneter Varianzschätzer für den ML-Schätzer ist dann \(\left(H(\beta)|_{\hat{\beta}_0,\hat{\beta}_1}\right)^{-1}\)

Ein halbgarer Ansatz für \(\boldsymbol{H}\):

\[\begin{aligned} \boldsymbol{H}(\beta) &= - \frac{\partial^2 l(\beta)}{\partial\beta\partial\beta^T} = \sum_{i=1}^n \boldsymbol{x}_i^T\boldsymbol{x}_i \pi_i (1 - \pi_i) \\ &= \sum_{i=1}^n \boldsymbol{x}_i^T\boldsymbol{x}_i \frac{\exp(\beta_0 + x \beta_1)}{1 + \exp(\beta_0 + x \beta_1)} \left(1 - \frac{\exp(\beta_0 + x \beta_1)}{1 + \exp(\beta_0 + x \beta_1)}\right) \\ &= \sum_{i=1}^n \boldsymbol{x}_i^T\boldsymbol{x}_i \frac{\exp(\beta_0 + x \beta_1)}{(1 + \exp(\beta_0 + x \beta_1))^2} \\ &= ?? \end{aligned}\]

Allerdings bekommen wir’s nicht wirklich vereinfacht, geschweige denn invertiert.