Aufgabe 1

Zeigen Sie die im Foliensatz “Kategoriale Regression” getätigte Aussage, dass die parametrische Familie der Multinomialverteilungen eine mehrdimensionale Exponentialfamilie ist. Bestimmen Sie über diese Eigenschaft den Erwartungswert und die Kovarianzmatrix einer multinomialverteilten Zufallsvariable

\[\begin{aligned} Y & \sim MN(m,\boldsymbol{\pi})\\ f(\boldsymbol{Y}\vert\boldsymbol{\pi}) &= \frac{m!}{\prod^q_{i=1}y_i!\cdot y_c!} \cdot \prod^q_{i=1}\pi_i^{y_i}\cdot \pi_c^{y_c}\\[1.5em] \log(f(\boldsymbol{Y}\vert\boldsymbol{\pi})) &= \log\underbrace{\left(\frac{m!}{\prod^q_{i=1}y_i! \cdot y_c!}\right)}_{c}+\sum^q_{i=1}\log\pi_i^{y_i}+\log\pi_c^{y_c}\\ &=\log c + \sum^q_{i=1}y_i\log\pi_i+ y_c\log\pi_c \end{aligned}\]

Mit \(\pi_c= 1-\sum^q_{i=1}\pi_i\) und \(y_c = m - \sum^q_{i=1}y_i\):

\[\begin{aligned} \log(f(\boldsymbol{Y}\vert\boldsymbol{\pi})) &=\log c + \sum^q_{i=1}y_i\log\pi_i+\left( m - \sum^q_{i=1}y_i\right)\log\pi_c\\ \log(f(\boldsymbol{Y}\vert\boldsymbol{\pi})) &=\log c + \sum^q_{i=1}y_i\log\pi_i - \sum^q_{i=1}y_i\log\pi_c\log\pi_c+ m\log\pi_c\\ \end{aligned}\]

Der mittlere Term lässt sich ausdrücken als

\[\begin{aligned} \sum^q_{i=1}y_i\log\pi_i - \sum^q_{i=1}y_i\log\pi_c &= y_1\log\pi_1-y_1\log\pi_c+\dots+y_q\log\pi_q-y_q\log\pi_c \\ &= \sum^q_{i=1}y_i(\log\pi_i-\log\pi_c) \\ &=\sum^q_{i=1}y_i\left(\log\frac{\pi_i}{\pi_c}\right) \\ &=\boldsymbol{Y\theta} \end{aligned}\]

\[\begin{aligned} \log(f(\boldsymbol{Y}\vert\boldsymbol{\pi})) &=\log c + \boldsymbol{Y\theta}+ m\log\pi_c\\ \end{aligned}\]

Für \(b(\theta)\):

\[\begin{aligned} \log\pi_c &= -\log\frac{1}{\pi_c}\\ &=-\log\frac{\pi_c+1-\pi_c}{\pi_c}\\ &=-\log\left(1+\frac{1-\pi_c}{\pi_c}\right)\\ &=-\log\left(1+\frac{\sum^q_{s=1}\pi_s}{\pi_c}\right)\\ &=-\log\left(1+\sum^q_{s=1}\frac{\pi_s}{\pi_c}\right)\\ &=-\log\left(1+\sum^q_{s=1}e^{\log\left(\frac{\pi_s}{\pi_c}\right)}\right)\\ &=-\log\left(1+\sum^q_{s=1}e^{\theta_s}\right)\\ &=-b(\boldsymbol\theta) \end{aligned}\]

Es ergibt sich also für \(f\) die Form einer Exponentialfamilie

\[\begin{aligned} f(\boldsymbol{Y}\vert\boldsymbol{\theta}, \phi, \omega) = \exp\left(\frac{\boldsymbol{Y\theta}-mb(\boldsymbol{\theta})}{\phi}\omega-c(\boldsymbol{Y},\phi,\omega)\right) \end{aligned}\]

mit oben definierten \(\theta, b, c\) und \(\phi = \omega = 1\).

Erwartungswert von \(\boldsymbol{Y}\):

\[\begin{aligned} E(Y_i)&=m\frac{\partial b(\theta)}{\partial \theta_i}\\ &=m\frac{e^{\theta_i}}{1+\sum_{s=1}^q e^{\theta_s}}\\ &=m\frac{\pi_i/\pi_c}{1+(1-\pi_c)/\pi_c}\\ &=m\frac{\pi_i/\pi_c}{\pi_c/\pi_c+(1-\pi_c)/\pi_c}\\ &=m\frac{\pi_i/\pi_c}{1/\pi_c}\\ &=m\pi_i\\ \Rightarrow E(\boldsymbol{Y})&=m\frac{\partial b(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}}=m\boldsymbol{\pi} \end{aligned}\]

Kovarianzmatrix:

\[\begin{aligned} \text{Cov}(\boldsymbol{Y})&=m\frac{\partial^2 b(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^T}=m\left(\frac{\partial^2 b(\boldsymbol{\theta})}{\partial \theta_r\partial \theta_s}\right)_{r,s=1}^q\\ &=m\frac{0\cdot(1+\sum_{i=1}^q e^{\theta_i})-e^{\theta_r}e^{\theta_s}}{\left(\sum_{i=1}^q e^{1+\theta_i}\right)^2}\\ &=-m\pi_r\pi_s \end{aligned}\]

Für die Varianz und damit die Hauptdiagonale der Kovarianzmatrix setze \(r=s\):

\[\begin{aligned} m\frac{\partial^2 b(\boldsymbol{\theta})}{\partial \theta_r\partial \theta_r} &=m\frac{1\cdot e^{\theta_r}-e^{\theta_r}e^{\theta_r}}{\left(\sum_{i=1}^q e^{1+\theta_i}\right)^2}\\ &=m\pi_r(1-\pi_r) \end{aligned}\]

Aufgabe 2

Der R-Datensatz pneumo aus der library faraway enthält Angaben über die Erkrankung an Pneumokoniose (Pneumonoultramicroscopicsilicovolcanoconiosis, Staublunge) von Bergarbeitern. Die Schwere der Pneumokoniose wurde hierbei in drei Kategorien (normal, mild und severe) eingeteilt. Es ist jeweils die Anzahl an Bergleuten in der jeweiligen Kategorie angegeben, sowie die Anzahl der Jahre, die die jeweiligen Bergleute unter Tage gearbeitet haben.

data("pneumo", package = "faraway")

pneumo <- pneumo %>%
  mutate(
    status = factor(status), # Rekodierung als factor
    status = relevel(status, ref = "normal"), # "normal" als Referenzkategorie (?),
    status_ord = factor(status, levels = c("normal", "mild", "severe"), ordered = TRUE)
  )

# Blick auf die Daten
pneumo %>%
  ggplot(aes(x = factor(year), y = Freq, fill = status)) +
  geom_col(alpha = .75, color = "black", position = "dodge") +
  # scale_x_continuous(breaks = seq(0, 70, 5)) +
  scale_y_continuous(breaks = scales::pretty_breaks()) +
  scale_fill_viridis_d(begin = .2, end = .8, option = "A") +
  labs(
    title = "Pneumonoultramicroscopicsilicovolcanoconiosis",
    subtitle = "Status nach Jahren",
    x = "Jahre", y = "Anzahl", fill = "Status"
  )

Wir sehen graphisch zumindest schonmal, dass “normal” am häufigsten vorkommt, und über die Zeit seltener wird, während zuerst “mild” und dann “severe” im Zeitverlauf an Häufigkeit zunehmen. Der Unterschied zwischen “mild” und “severe” ist dabei relativ klein – darauf nehmen wir später nochmal Bezug.

a)

Behandeln Sie den Pneumokoniose-Status als nominal und passen Sie ein mehrkategorielles Logit-Modell mit der Kovariable years an. Berechnen Sie die prognostizierte Verteilung der Status-Variable für einen Bergarbeiter, der 25 Jahre gearbeitet hat.

# Model fit
mod_nom <- nnet::multinom(status ~ year, weights = Freq, data = pneumo, trace = FALSE)

# Gleiches Modell, aber ohne weights, dafür mit data-reshaping - gleiche Ergebnisse
# faraway::pneumo %>%
#   tidyr::spread(status, Freq, fill = 0) %>%
#   nnet::multinom(cbind(normal, mild, severe) ~ year, data = .)

# Model coefficients
table_nom <- mod_nom %>%
  broom::tidy(exponentiate = FALSE) %>%
  mutate(OR = exp(estimate), p.value = ifelse(p.value < 0.05, "< 0.05", p.value)) %>%
  select(y.level, term, estimate, std.error, OR, statistic, p.value) 

table_nom %>%
  kable(
    booktabs = TRUE, digits = 3,
    col.names = c("Status", "Term", "Koeffizient", "SE", "OR", "z", "p")
  ) %>%
  kable_styling(position = "center")
Status Term Koeffizient SE OR z p
mild (Intercept) -4.292 0.521 0.014 -8.231 < 0.05
mild year 0.084 0.015 1.087 5.469 < 0.05
severe (Intercept) -5.060 0.596 0.006 -8.484 < 0.05
severe year 0.109 0.016 1.115 6.636 < 0.05
# Predicted value: year = 25
predict(mod_nom, newdata = data.frame(year = 25))
## [1] normal
## Levels: normal mild severe

Wir erhalten Koeffizienten für die Kategorien “mild” und “severe”, jeweils mit Referenzkategorie “normal”. Die Intercepts ignorierend erhalten wir:

  • mild: Ein zusätzliches Jahr unter Tage (year) ist mit einem um den Faktor ~1.09 höheren Wahrscheinlichkeit assoziiert, einen leichten Fall von Staublunge zu entwickeln.
  • severe: Analog “mild” ist ist hier der Faktor ~1.12, mit dem sich die Wahrscheinlichkeit einen schwerwiegenden Fall Staublunge zu entwickeln mit einem zusätzlichen Jahr unter Tage erhöht.

Unsere Prognose für einen Bergarbeitr, der 25 Jahre unter Tage gearbeitet hat, wäre in diesem Modell “normal”, also (noch) kein Fall von Staublunge.

b)

Führen Sie die gleiche Analyse wie in a) durch. Behandeln Sie jedoch diesmal die Status-Variable als ordinal, indem Sie ein kumulatives Logit-Modell anpassen.

# Model fit
mod_ord <- MASS::polr(status_ord ~ year, weights = Freq, data = pneumo)

# Alternative:
# ordinal::clm(status_ord ~ year, weights = Freq, data = pneumo, link = "logit")

# Model coefficients
table_ord <- mod_ord %>%
  broom::tidy(exponentiate = FALSE) %>%
  mutate(OR = exp(estimate)) %>%
  select(coefficient_type, term, estimate, std.error, OR, statistic) 

table_ord %>%
  kable(
    booktabs = TRUE, digits = 3,
    col.names = c("Koeffizient-Art", "Term", "Koeffizient", "SE", "exp(Koef)","z")
  ) %>%
  kable_styling(position = "center")
Koeffizient-Art Term Koeffizient SE exp(Koef) z
coefficient year 0.096 0.012 1.101 8.034
zeta normal|mild 3.956 0.410 52.240 9.656
zeta mild|severe 4.869 0.441 130.197 11.038
# Predicted value: year = 25
predict(mod_ord, newdata = data.frame(year = 25))
## [1] normal
## Levels: normal mild severe

Wir betrachten hier Status als geordnet: \(\mathtt{normal} < \mathtt{mild} < \mathtt{severe}\).

Wir erhalten einen Koeffizienten für year, der uns sagt, dass sich das Risiko in die “nächst schlechtere” Kategorie Staublunge aufzusteigen pro zusätzlichem Jahr um den Faktor ~1.1 erhöht.

Die beiden “zeta” Koeffizienten sind dabei wohl die Schwellenwerte aus dem Modell, und sollten anzeigen, dass ab ~52 Jahren der “Übergang” von “normal” zu “mild” stattfindet. Für 25 Jahre prognostiziert das Modell “normal”, was sich damit decken würde. Der Übergang zu “severe” findet erst bei etwas unrealistischen ~130 Jahren statt, was vermutlich an der eingangs anhand des Plots veranschaulichten schwachen Differenzierung zwischen den beiden Krankheitskategorien liegt.

Aus Übung

\[\begin{aligned} P(Y_i \leq r) &= \frac{\exp(\theta_r + x_i \beta)}{1 + \exp(\theta_r + x_i \beta)} \\ &\Rightarrow \log \left(\frac{P(y_i \leq r |x_i)}{P(y_i > r |x_i)}\right) = \theta_r + x_i \beta \end{aligned}\]

\(\exp(0.096) = 1.108\) Log-Odds der kum. Wahrsk. steigt um 0.096.
Odds der kum. Wahrsk. steigt um ca 10%.

c)

Passen Sie nun ein kumulatives Extremwert-Modell an die Daten an und wiederholen Sie die Analyse.

# Model fit (link = cumulative log-log -> method = "cloglog")
mod_cloglog <- MASS::polr(status_ord ~ year, method = "cloglog", weights = Freq, data = pneumo)

# Alternative:
# ordinal::clm(status_ord ~ year, weights = Freq, data = pneumo, link = "cloglog")

# Model coefficients
table_cloglog <- mod_cloglog %>%
  broom::tidy(exponentiate = FALSE) %>%
  mutate(OR = exp(estimate)) %>%
  select(coefficient_type, term, estimate, std.error, OR, statistic) 

table_cloglog %>%
  kable(
    booktabs = TRUE, digits = 3,
    col.names = c("Koeffizient-Art", "Term", "Koeffizient", "SE", "exp(Koef)", "z")
  ) %>%
  kable_styling(position = "center")
Koeffizient-Art Term Koeffizient SE exp(Koef) z
coefficient year 0.050 0.006 1.052 8.362
zeta normal|mild 1.746 0.175 5.731 9.987
zeta mild|severe 2.220 0.188 9.212 11.794
# Predicted value: year = 25
predict(mod_cloglog, newdata = data.frame(year = 25))
## [1] normal
## Levels: normal mild severe

Ähnlich wie beim vorigen Modell, allerdings ist nun ~1.05 der Faktor mit dem sich die Wahrscheinlichkeit ändert, die Schwellenwerte zu überschreiten.

Aus Übung:

Wir modellieren die negativ-logarithmierte Übertretungswahrscheinlichkeit.

-> Schwer zu interpretieren :(

d)

Vergleichen und interpretieren Sie die drei Modelle.

bind_rows(
  table_nom %>% 
    mutate(
      Modell = "Nominal", 
      coefficient_type = y.level
    ) %>%
    select(-y.level, -p.value),
  table_ord %>% mutate(Modell = "Ordinal"),
  table_cloglog %>% mutate(Modell = "Extremwert (cloglog)")
) %>%
  select(Modell, coefficient_type, everything()) %>%
  mutate_if(is.numeric, round, 3) %>%
  kable(
    booktabs = TRUE,
    col.names = c("Modell", "Koeffizient-Typ", "Term", "Koeffizient", "SE", "exp(Koef)", "z")
  ) %>%
  kable_styling(position = "center") %>%
  column_spec(column = 1, bold = TRUE) %>%
  collapse_rows(columns = 1)
Modell Koeffizient-Typ Term Koeffizient SE exp(Koef) z
Nominal mild (Intercept) -4.292 0.521 0.014 -8.231
mild year 0.084 0.015 1.087 5.469
severe (Intercept) -5.060 0.596 0.006 -8.484
severe year 0.109 0.016 1.115 6.636
Ordinal coefficient year 0.096 0.012 1.101 8.034
zeta normal|mild 3.956 0.410 52.240 9.656
zeta mild|severe 4.869 0.441 130.197 11.038
Extremwert (cloglog) coefficient year 0.050 0.006 1.052 8.362
zeta normal|mild 1.746 0.175 5.731 9.987
zeta mild|severe 2.220 0.188 9.212 11.794

Das nominale Modell sagt uns, dass das Risiko einen mild/severe Fall Staublunge zu entwickeln mit der Zeit unter Tage steigt, uns zwar für beide Kategorien getrennt, ohne Berücksichtigung der Ordnungsbeziehung zwischen den beiden Schweregraden.

Die kumulativen Modelle berücksichtigen diese, und sparen uns einen Parameter. Das ordinale Modell sagt uns, dass das Risiko in die nächste Kategorie überzugehen mit einem Jahr um den Faktor ~1.1 steigt, wobei uns das Extremwertmodell sagt, dass sich die kategorieunabhängige Schwellenwertüberschreitungswahrscheinlichkeit mit einem zusätzlichen Jahr um den Faktor ~1.05 ändert.