Aufgabe 1

Zwischen Januar 1972 und April 1974 wurde im französischen Département Ille-et-Vilaine (Bretagne) eine Fall-Kontroll-Studie zu Risikofaktoren des Ösophagus-Krebses durchgeführt. 200 männliche Fälle und 775 Kontrollen konnten letztendlich zu ihren Ernährungsgewohnheiten befragt werden. Die Kontrollen wurden zufällig aus Wählerverzeichnissen gezogen. Daten dieser Studie sollen hier untersucht werden. Für eine detaillierte Studienbeschreibung sei auf Breslow and Day (1980) und Tuyns et al. (1977) verwiesen.

Der Datensatz zu dieser Studie (ille_et_vilaine.csv) beinhaltet aggregierte Informationen zu Alter, Alkohol- und Nikotinkonsum der Studienteilnehmer.

  • Alter (Variable alt) liegt in 6 Kategorien vor (1=25-34 Jahre, 2=35-44 Jahre, 3=45-54 Jahre, 4=55-64 Jahre, 5=65-74 Jahre, 6= 75+ Jahre),

  • Alkoholkonsum (Variable alk) in 4 Kategorien (1=0-39 g/Tag, 2=40-79 g/Tag, 3=80-119 g/Tag, 4=120+ g/Tag) und

  • Nikotinkonsum (Variable tab) in 4 Kategorien (1=0-9 g/Tag, 2=10-19 g/Tag, 3=20-29 g/Tag, 4=30+ g/Tag).

Zu jeder möglichen Variablenkombination liegen die Anzahl an Beobachtungen (Variable n) sowie die Anzahl an Fällen (Variable n1) vor.

Versuchen Sie ein möglichst passendes logistisches Modell für die Daten zu finden:

Wir kodieren die UVs als factor, um sie nominalskaliert zu betrachten:

lille <- read_csv("data/ille_et_vilaine.csv", col_types = cols()) %>%
  mutate(alt = factor(alt),
         alk = factor(alk),
         tab = factor(tab))

a)

Passen Sie fünf Modelle mit den folgenden Einflussvariablen an die Daten an:

  • Modell 1: Alter
  • Modell 2: Alter + Nikotinkonsum
  • Modell 3: Alter + Alkoholkonsum
  • Modell 4: Alter + Nikotinkonsum + Alkoholkonsum
  • Modell 5: Alter + Nikotinkonsum + Alkoholkonsum + Interaktionen zwischen Nikotin und Alkoholkonsum
# n1 = Fälle ---  n - n1 = Kontrollen
m1 <- glm(cbind(n1, n - n1) ~ alt,                       family = binomial, data = lille)
m2 <- glm(cbind(n1, n - n1) ~ alt + tab,                 family = binomial, data = lille)
m3 <- glm(cbind(n1, n - n1) ~ alt + alk,                 family = binomial, data = lille)
m4 <- glm(cbind(n1, n - n1) ~ alt + tab + alk,           family = binomial, data = lille)
m5 <- glm(cbind(n1, n - n1) ~ alt + tab + alk + alk:tab, family = binomial, data = lille)

b)

Beurteilen Sie anhand der Deviance, der Pearson-Statistik und des AIC, welches der Modelle am geeignetsten ist.

Hinweis 1: Es seien zwei hierarchisch geordnete logistische Regressionmodelle M1 (mit \(p\) Parametern) und M2 (mit \(q\) Parametern), \(p<q\), gegeben. Unter der Nullhypothese, dass die zusätzlichen \(q-p\) Parameter von M2 alle Null sind, ist die Differenz der Deviance von M1 und M2 asymptotisch \(\chi^2\) -verteilt mit \(q-p\) Freiheitsgeraden.

# Extract model stats from each model, combine to table
list(m1, m2, m3, m4, m5) %>%
  purrr::map_df(~{
    tibble::tibble(
      Modell = as.character(formula(.x))[[3]], # UVs
      df = length(.x$coefficients), # Freiheitsgrade (inkl. intercept)
      Deviance = deviance(.x),
      Pearson = sum(residuals(.x, type = "pearson")^2),
      AIC = AIC(.x)
    )
  }) %>%
  kable(booktabs = TRUE, digits = 3) %>%
  kable_styling(position = "center")
Modell df Deviance Pearson AIC
alt 6 246.909 363.909 373.964
alt + tab 9 210.270 348.942 343.325
alt + alk 9 105.881 117.501 238.936
alt + tab + alk 12 82.337 86.557 221.392
alt + tab + alk + alk:tab 21 76.886 78.125 233.941

Nach Pearson-Statistik

# Berechnung der Pearson-Statistik: Summe der quadrierten Pearson-Residuen
pearson <- function(mod) {
  sum(residuals(mod, type = "pearson")^2)
}

#' Funktion zum Modellvergleich
#' @param x,y glm-Objekte, wobei x das "kleinere" Modell ist
#' @param compare_method Funktion zu Berechnung der Teststatistik (ohne "")
glm_compare <- function(x, y, compare_method = c(pearson, deviance)) {
  xname <- deparse(substitute(x))
  yname <- deparse(substitute(y))
  p1 <- compare_method(x)
  p2 <- compare_method(y)
  # Difference in statistics
  stat <- p1 - p2
  # Difference in dof, number of coefs minus intercept
  df_diff <- length(y$coefficients) - length(x$coefficients)
  # p value based on chi^2 assumption
  p <- pchisq(stat, df = df_diff, lower.tail = FALSE)
  p <- ifelse(p < 0.05, "< 0.05", as.character(round(p, 3)))
  # Return compactly
  tibble(
    comp = paste(xname, yname, sep = " vs. "), 
    stat = round(stat, 2), df = df_diff, p = p
  )
  # glue::glue("Differenz: {round(stat, 3)} (df = {df_diff}), p = {p}")
}

bind_rows(
  glm_compare(m1, m2, pearson),
  glm_compare(m1, m3, pearson),
  glm_compare(m2, m4, pearson),
  glm_compare(m3, m4, pearson),
  glm_compare(m4, m5, pearson)
) %>%
  kable(
    booktabs = TRUE, col.names = c("Vergleich", "$\\Delta \\chi^2$", "df", "p"),
    escape = FALSE
  ) %>%
  kable_styling(protect_latex = TRUE, position = "center")
Vergleich \(\Delta \chi^2\) df p
m1 vs. m2 14.97 3 < 0.05
m1 vs. m3 246.41 3 < 0.05
m2 vs. m4 262.38 3 < 0.05
m3 vs. m4 30.94 3 < 0.05
m4 vs. m5 8.43 9 0.491

Die hinzunahme weiterer Parameter führt zu einer signifikanten Reduktion der Pearson-Statistik bis zum Vergleich von Modell 3 und Modell 4, das heißt wir “stoppen” bei Modell 3, weil der zusätzliche Term (Interkationseffekt aus Alkohol + Nikotin) die Deviance nicht mehr signifikant reduziert.

Nach Deviance

bind_rows(
  glm_compare(m1, m2, deviance),
  glm_compare(m1, m3, deviance),
  glm_compare(m2, m4, deviance),
  glm_compare(m3, m4, deviance),
  glm_compare(m4, m5, deviance)
) %>%
  kable(
    booktabs = TRUE, col.names = c("Vergleich", "$\\Delta$ Deviance", "df", "p"),
    escape = FALSE
  ) %>%
  kable_styling(position = "center", protect_latex = TRUE)
Vergleich \(\Delta\) Deviance df p
m1 vs. m2 36.64 3 < 0.05
m1 vs. m3 141.03 3 < 0.05
m2 vs. m4 127.93 3 < 0.05
m3 vs. m4 23.54 3 < 0.05
m4 vs. m5 5.45 9 0.793

Wir kommen zum selben Schluss wie im vorigen Teil mit der Pearson-Statistik

Nach AIC:

Modell 4 (Alter, Nikotin, Tabak) hat das kleinste AIC und ist demnach zu bevorzugen.

Aufgabe 2

Die Dichte der inversen Gaußverteilung (\(Y \sim IG(\mu, \sigma^2)\)) ist gegeben durch

\[\begin{aligned} f(y) = \begin{cases} \left(\frac{\sigma^2}{2\pi y^3}\right) e^{-\frac{\sigma^2(y-\mu)^2}{2\mu^2 y}}, & y > 0 \\ 0, & y \leq 0 \end{cases} \end{aligned}\]

Zeigen Sie, dass die inversen Gaußverteilungen eine Exponentialfamilie bilden und bestimmen Sie Erwartungswert und Varianz der inversen Gaußverteilung.

\[\begin{aligned} f(y\ |\ \mu, \sigma^2) &= \left(\frac{\sigma^2}{2\pi y^3}\right)^{\frac{1}{2}} e^{-\frac{\sigma^2(y-\mu)^2}{2\mu^2 y}} && \big| \log \\ &= \log\left(\left(\frac{\sigma^2}{2\pi y^3}\right)^{\frac{1}{2}}\right) - \frac{\sigma^2(y-\mu)^2}{2\mu^2 y} \\ &= - \frac{\sigma^2(y-\mu)^2}{2\mu^2 y} + \frac{1}{2}\log\left(\frac{\sigma^2}{2\pi y^3}\right) \\ &= - \frac{\sigma^2(y^2 - 2y\mu + \mu^2)}{2\mu^2 y} + \frac{1}{2}\log\left(\frac{\sigma^2}{2\pi y^3}\right) \\ &= - \frac{\sigma^2 y^2 - \sigma^2 2y\mu + \sigma^2 \mu^2}{2\mu^2 y} + \frac{1}{2}\log\left(\frac{\sigma^2}{2\pi y^3}\right) \\ &= - \frac{\sigma^2 y}{2\mu^2} + \frac{\sigma^2}{\mu} - \frac{\sigma^2}{2y} + \frac{1}{2}\log\left(\frac{\sigma^2}{2\pi y^3}\right) \\ &= - \frac{\sigma^2}{2\mu^2} y + \frac{\sigma^2}{\mu} + \underbrace{ \frac{1}{2}\log\left(\frac{\sigma^2}{2\pi y^3}\right) - \frac{\sigma^2}{2y}}_ {:=\ c(y,\ \phi,\ \omega\ = 1)} \\ &= - \frac{\theta \sigma^2}{2} y + \theta^{\frac{1}{2}}\sigma^2 + c(y, \phi, \omega) && \big|\ \text{sei } \theta = \frac{1}{\mu^2} \Rightarrow \mu = \theta^{-\frac{1}{2}}\\ &= \frac{\theta y}{\phi} + \frac{-2\theta^{\frac{1}{2}}}{\phi} + c(y, \phi, \omega) && \big|\ \text{sei } \phi = \frac{-2}{\sigma^2} \Rightarrow \sigma^2 = \frac{-2}{\phi}\\ &= \frac{\theta y - 2\theta^{\frac{1}{2}}}{\phi} + c(y, \phi, \omega) && \big|\ \text{sei } b(\theta) = 2\theta^{\frac{1}{2}} \\ f(y\ |\ \theta, \phi, \omega) &= \exp\left\{\frac{\theta y - b(\theta)}{\phi} + c(y, \phi, \omega)\right\} \\&&\Box \end{aligned}\]

Erwartungswert und Varianz, für letztere ausgehend von \(\omega = 1\):

\[\begin{aligned} b(\theta) &= 2\theta^{\frac{1}{2}} \\ \mathbb{E}(Y) &= b'(\theta) = \theta^{-\frac{1}{2}} = \left(\frac{1}{\mu^2}\right)^{-\frac{1}{2}} = \mu \\ \mathbb{V}(Y) &= \phi \cdot b''(\theta) = \frac{-2}{\sigma^2} \cdot \frac{-1}{2} \cdot \theta^{-\frac{3}{2}} = \frac{\mu^3}{\sigma^2} \end{aligned}\]