Lösungen

Vergleiche zuerst mit deinem eigenen Versuch. Es gibt oft mehrere richtige Wege.

Aufgabe 1 – Einfache logistische Regression

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

df = pd.read_csv("https://schradern.github.io/data-science-coach/data/kohorte.csv")
result = smf.logit("verstorben_30d ~ sofa_score", data=df).fit(disp=False)
ci = result.conf_int()
print("OR sofa_score:   ", round(np.exp(result.params["sofa_score"]), 3))
print("95% CI lower:    ", round(np.exp(ci.loc["sofa_score", 0]), 3))
print("95% CI upper:    ", round(np.exp(ci.loc["sofa_score", 1]), 3))

suppressPackageStartupMessages(library(tidyverse))

df <- read_csv("https://schradern.github.io/data-science-coach/data/kohorte.csv", show_col_types = FALSE)
m  <- glm(verstorben_30d ~ sofa_score, data = df, family = binomial)
round(exp(coef(m)), 3)
round(exp(suppressMessages(confint(m))), 3)

Interpretation: Eine OR von ~1,60 (95-%-KI ≈ [1,41; 1,82]) bedeutet: jeder zusätzliche SOFA-Punkt geht mit ~60 % höheren Chancen (Odds) auf 30-Tage-Mortalität einher. Das ist ein unadjustierter Schätzer; alle anderen Risikofaktoren, die mit dem SOFA-Score korrelieren (Alter, Aufnahmegrund), sind noch nicht herausgerechnet.

Der wahre unadjustierte OR beträgt 1,477, der wahre adjustierte OR 1,387 (beide aus lib/ground_truth.py). Das einfache Modell zielt also von vornherein auf eine größere Zahl, und die Stichprobe streut zusätzlich darum.

Aufgabe 2 – Multivariables Modell und Wahrheitsvergleich

from lib.ground_truth import true_odds_ratios_for

df["sepsis"] = (df["aufnahmegrund"] == "Sepsis").astype(int)
df["age_z"]  = df["alter"] - 64
formula = "verstorben_30d ~ age_z + sofa_score + crp_mg_l + diabetes + sepsis"
result2 = smf.logit(formula, data=df).fit(disp=False)

wahr = true_odds_ratios_for(("age_centred", "sofa_score", "crp_mg_l", "diabetes", "sepsis"))
print(np.exp(result2.params).round(3))
print(wahr)

df <- df |> mutate(sepsis = as.integer(aufnahmegrund == "Sepsis"),
                   age_z  = alter - 64)
m2 <- glm(verstorben_30d ~ age_z + sofa_score + crp_mg_l + diabetes + sepsis,
          data = df, family = binomial)
round(exp(coef(m2)), 3)

Beobachtung: Das 95-%-KI überdeckt die wahre OR bei allen fünf Prädiktoren. Es liegt also keine Verzerrung vor, nur Stichprobenvarianz.

Relativ am stärksten weicht diabetes ab (1,11 vs. 1,70). Der Grund ist Präzision, nicht Bias: bei 78 Ereignissen und ~25 % Diabetes-Prävalenz reicht das KI von 0,58 bis 2,14 — es enthält sowohl „schützend" als auch „verdoppelt das Risiko". crp_mg_l sieht relativ ebenfalls schlecht aus, weil sein wahrer Effekt winzig ist; jeder absolute Fehler wirkt dort relativ groß.

Der wichtigste Punkt der Aufgabe: Die wahren Werte sind 1,053 / 1,384 / 1,004 / 1,703 / 2,036 — nicht exp(0,045) / exp(0,28) / …. Letztere sind Hazard Ratios (1,046 / 1,320 / 1,004 / 1,549 / 1,835) und liegen systematisch näher an 1. Bei einer Ereignisrate von 15,6 % ist die Näherung „OR ≈ HR" nicht mehr tragfähig. Wer sie trotzdem benutzt, produziert eine scheinbare Verzerrung von 5–12 % in dieselbe Richtung für jeden Prädiktor.

Aufgabe 3 – Confounder und Mediator unterscheiden

for exposure in ["sepsis", "diabetes"]:
    for label, extra in [("crude", ""), ("+ alter", " + age_z"),
                         ("+ alter + sofa", " + age_z + sofa_score")]:
        m = smf.logit(f"verstorben_30d ~ {exposure}{extra}", data=df).fit(disp=False)
        ci = np.exp(m.conf_int().loc[exposure])
        print(f"{exposure:9s} {label:16s} OR={np.exp(m.params[exposure]):.3f}"
              f"  [{ci[0]:.2f}, {ci[1]:.2f}]")

for (exposure in c("sepsis", "diabetes")) {
  for (extra in c("", " + age_z", " + age_z + sofa_score")) {
    m <- glm(as.formula(paste0("verstorben_30d ~ ", exposure, extra)),
             data = df, family = binomial)
    cat(sprintf("%-9s %-20s OR=%.3f\n", exposure, extra, exp(coef(m)[[exposure]])))
  }
}

Sepsis — hier ist der Bias sichtbar

alter ist ein Confounder. Ältere Patient:innen sind kränker und sterben unabhängig davon häufiger. Für ihn muss adjustiert werden — begründet aus der DAG, nicht daraus, dass sich die Zahl ändert.

sofa_score ist ein Mediator. Laut DAG (data/README.md) liegt er auf dem Pfad sepsis → sofa_score → Tod: eine Sepsis erhöht den SOFA-Score um ~3 Punkte, jeder Punkt erhöht die Hazard. Adjustiert man dafür, kontrolliert man genau den Teil des Effekts weg, nach dem man gefragt hat. Die OR fällt von 5,75 auf 1,92, und die Konfidenzintervalle überlappen kaum — der Overadjustment-Bias ist hier klar belegt, nicht nur behauptet.

Was du der Klinikleitung berichtest: den totalen Effekt (OR ≈ 5,8; wahr 4,68), also das Modell ohne sofa_score. Die Frage „Wie viel gefährlicher ist eine Sepsis?" schließt den Weg über den Organversagensscore ausdrücklich mit ein. Für eine Mediationsanalyse („Wie viel des Effekts läuft über SOFA?") ist umgekehrt das SOFA-adjustierte Modell das richtige Werkzeug (→ Modul 15).

Diabetes — dasselbe Muster, aber nicht belegbar

Die Punktschätzer fallen genauso — 2,07 → 1,59 → 0,99. Aber jedes dieser Konfidenzintervalle enthält sowohl den wahren totalen als auch den wahren direkten Effekt. Bei N = 500 und 78 Ereignissen lässt sich der Overadjustment-Bias für Diabetes herleiten (aus der DAG), aber nicht zeigen (aus den Daten). Der Sprung auf 0,99 ist zum großen Teil Rauschen.

Aufgabe 4 – Lineare Regression und Residuenplot

import matplotlib.pyplot as plt
import statsmodels.formula.api as smf

m_ols = smf.ols("verweildauer_tage ~ sofa_score + age_z", data=df).fit()
fitted    = m_ols.fittedvalues
residuals = m_ols.resid

plt.figure(figsize=(7, 4))
plt.scatter(fitted, residuals, alpha=0.4, s=20)
plt.axhline(0, color="red", linewidth=1)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.title("Residual plot — OLS for Verweildauer")
plt.tight_layout()
plt.savefig("residualplot.png", dpi=120)
print(f"R²: {m_ols.rsquared:.4f}")

m_ols <- lm(verweildauer_tage ~ sofa_score + age_z, data = df)
plot(fitted(m_ols), resid(m_ols),
     main = "Residualplot — OLS Verweildauer",
     xlab = "Fitted values", ylab = "Residuals")
abline(h = 0, col = "red")

Befund: Der Residuenplot zeigt typischerweise eine trichterförmige Streuung (Heteroskedastizität) und rechtschiefe Residuen, Zeichen, dass OLS die Verweildauer-Verteilung schlecht abbildet. Ein Gamma-GLM modelliert den strikt positiven, rechtsschiefen Outcome direkt:

import statsmodels.api as sm          # die Familien liegen in statsmodels.api,
import statsmodels.formula.api as smf  # nicht in statsmodels.formula.api

m_gamma = smf.glm("verweildauer_tage ~ sofa_score + age_z", data=df,
                  family=sm.families.Gamma(link=sm.families.links.Log())).fit()
print(m_gamma.summary())

Aufgabe 5 – Cox-Regression

from lifelines import CoxPHFitter

# fu_zeit_tage/status are the real time-to-event pair — NOT verweildauer_tage
# (that column is descriptive length of stay, see data/README.md).
cox_df = df[["fu_zeit_tage", "status", "sofa_score", "sepsis"]].copy()
cph = CoxPHFitter()
cph.fit(cox_df, duration_col="fu_zeit_tage", event_col="status")
cols = ["exp(coef)", "exp(coef) lower 95%", "exp(coef) upper 95%", "p"]
print(cph.summary[cols])

suppressPackageStartupMessages(library(survival))
m_cox <- coxph(Surv(fu_zeit_tage, status) ~ sofa_score + sepsis,
               data = df)
summary(m_cox)$conf.int

Interpretation: Sepsis: HR ≈ 1,51 (95-%-KI ≈ [0,88; 2,59], p ≈ 0,14); SOFA-Score: HR ≈ 1,43 (95-%-KI ≈ [1,27; 1,61], p < 0,0001). Eine HR von ~1,5 für Sepsis bedeutet: septische Patient:innen haben zu jedem Zeitpunkt ein ~50 % höheres momentanes Sterberisiko als Patient:innen mit anderem Aufnahmegrund, unter der Voraussetzung, dass die proportionale-Hazards-Annahme gilt (d. h. dieses Verhältnis bleibt über die gesamte Beobachtungszeit konstant). Bei diesem einfachen Modell (nur SOFA + Sepsis, N = 500, 78 Ereignisse) ist der Sepsis-Effekt für sich genommen nicht signifikant, das KI ist breit und schließt 1 ein, die Punktschätzung passt aber zur dokumentierten Ground Truth (Sepsis-HR ≈ 1,56, siehe data/README.md). Prüfung der PH-Annahme mit cox.zph() (R) bzw. cph.check_assumptions() (lifelines).

Bonus – EPV und Modellkomplexität

n_events = int(df["verstorben_30d"].sum())   # 78 at N=500 (15.6 % mortality)
n_pred   = 8
epv      = n_events / n_pred
print(f"EPV: {epv:.2f}")
# EPV just under 10 → model with 8 predictors is borderline for this cohort

Diskussion: Mit 78 Ereignissen und 8 Prädiktoren liegt EPV = 9,75, knapp unterhalb der klassischen Faustregel (≥ 10). Mit den 6 Prädiktoren des Hauptmodells in diesem Modul liegt EPV dagegen bei ≈ 13 (solide vertretbar), das Bonus-Szenario mit 8 Prädiktoren zeigt also konkret, wie schnell zusätzliche Variablen die Faustregel unterschreiten. Die Schätzungen sind bei EPV ≈ 9,75 noch möglich, aber mit breiteren Konfidenzintervallen und erhöhtem Overfitting-Risiko verbunden als beim Hauptmodell.

Die Faustregel ist allerdings grob. Riley et al. (2019, Statistics in Medicine) zeigen, dass die benötigte Stichprobengröße auch von der erwarteten Effektgröße, dem R² des Gesamtmodells und der Ereignisrate abhängt. Dein pmsampsize-Paket (R) berechnet die tatsächlich benötigte N für ein gegebenes Modell, eine präzisere Alternative zur EPV-Faustregel.