Lösungen

import pandas as pd
import numpy as np
from scipy import stats

cohort = pd.read_csv("https://schradern.github.io/data-science-coach/data/kohorte.csv")
labs   = pd.read_csv("https://schradern.github.io/data-science-coach/data/labor.csv")
df     = cohort.merge(labs, on="patient_id", how="left")

with_diabetes    = df.loc[df["diabetes"] == 1, "sofa_score"].dropna()
without_diabetes = df.loc[df["diabetes"] == 0, "sofa_score"].dropna()

t_stat, p = stats.ttest_ind(with_diabetes, without_diabetes, equal_var=False)

# Cohen's d (pooled SD)
pooled_sd = np.sqrt(
    ((len(with_diabetes) - 1) * with_diabetes.std()**2 +
     (len(without_diabetes) - 1) * without_diabetes.std()**2) /
    (len(with_diabetes) + len(without_diabetes) - 2)
)
d = (with_diabetes.mean() - without_diabetes.mean()) / pooled_sd

# 95-%-CI of the difference (Welch approximation)
se_diff = np.sqrt(with_diabetes.var(ddof=1) / len(with_diabetes) +
                  without_diabetes.var(ddof=1) / len(without_diabetes))
df_welch = (with_diabetes.var(ddof=1) / len(with_diabetes) +
            without_diabetes.var(ddof=1) / len(without_diabetes))**2 / (
    (with_diabetes.var(ddof=1) / len(with_diabetes))**2 / (len(with_diabetes) - 1) +
    (without_diabetes.var(ddof=1) / len(without_diabetes))**2 / (len(without_diabetes) - 1))
t_crit   = stats.t.ppf(0.975, df=df_welch)
diff     = with_diabetes.mean() - without_diabetes.mean()
ci_lo, ci_hi = diff - t_crit * se_diff, diff + t_crit * se_diff

print(f"Diabetes:         n={len(with_diabetes)}, SOFA={with_diabetes.mean():.2f}")
print(f"kein Diabetes:    n={len(without_diabetes)}, SOFA={without_diabetes.mean():.2f}")
print(f"Welch-t={t_stat:.3f}, p={p:.4f}, Cohen's d={d:.3f}")
print(f"Differenz={diff:.2f}, 95%-KI=[{ci_lo:.2f}, {ci_hi:.2f}]")

suppressPackageStartupMessages(library(tidyverse))

cohort <- read_csv("https://schradern.github.io/data-science-coach/data/kohorte.csv", show_col_types = FALSE)
labs   <- read_csv("https://schradern.github.io/data-science-coach/data/labor.csv", show_col_types = FALSE)
df     <- left_join(cohort, labs, by = "patient_id")

with_d    <- df |> filter(diabetes == 1) |> pull(sofa_score) |> na.omit()
without_d <- df |> filter(diabetes == 0) |> pull(sofa_score) |> na.omit()

t_res <- t.test(with_d, without_d, var.equal = FALSE)
print(t_res)

# Cohen's d
pooled_sd <- sqrt(((length(with_d)-1)*var(with_d) + (length(without_d)-1)*var(without_d)) /
                  (length(with_d) + length(without_d) - 2))
d <- (mean(with_d) - mean(without_d)) / pooled_sd
cat(sprintf("Cohen's d = %.3f\n", d))

laktat = df["laktat_mmol_l"].dropna()
w_stat, p_sw = stats.shapiro(laktat)
print(f"Shapiro-Wilk: W={w_stat:.4f}, p={p_sw:.2e}")
print(f"n={len(laktat)}, Schiefe={laktat.skew():.2f}")
print(f"Normalverteilung abgelehnt (p<0.05): {p_sw < 0.05}")

laktat <- df$laktat_mmol_l |> na.omit()
sw_res <- shapiro.test(laktat)
print(sw_res)
cat(sprintf("Schiefe (approx.) = %.2f\n",
            (mean(laktat) - median(laktat)) / sd(laktat)))

import pandas as pd

# Binarize smoking status: active vs. former/never
df2 = df.copy()
df2["smoker_active"] = (df2["raucherstatus"] == "aktiv").astype(int)

ct = pd.crosstab(df2["smoker_active"], df2["verstorben_30d"],
                 rownames=["smoker_active"], colnames=["verstorben_30d"])
print(ct)

chi2, p_chi2, dof, _ = stats.chi2_contingency(ct, correction=False)
cramers_v = np.sqrt(chi2 / (ct.values.sum() * (min(ct.shape) - 1)))

print(f"χ²={chi2:.3f}, df={dof}, p={p_chi2:.4f}, Cramér's V={cramers_v:.3f}")

df2 <- df |> mutate(smoker_active = if_else(raucherstatus == "aktiv", 1L, 0L))
ct  <- table(smoker_active = df2$smoker_active, verstorben_30d = df2$verstorben_30d)
print(ct)

chi_res   <- chisq.test(ct, correct = FALSE)
print(chi_res)
cramers_v <- sqrt(chi_res$statistic / (sum(ct) * (min(dim(ct)) - 1)))
cat(sprintf("Cramér's V = %.3f\n", cramers_v))

groups = ["Sepsis", "Pneumonie", "Herzinsuffizienz"]

for group in groups:
    g  = df.loc[df["aufnahmegrund"] == group, "laktat_mmol_l"].dropna()
    n  = len(g)
    se = g.std() / np.sqrt(n)
    ci = stats.t.interval(0.95, df=n - 1, loc=g.mean(), scale=se)
    print(f"{group}: n={n}, Mittel={g.mean():.2f}, 95%-KI=[{ci[0]:.2f}, {ci[1]:.2f}]")

groups <- c("Sepsis", "Pneumonie", "Herzinsuffizienz")
for (gr in groups) {
  g  <- df |> filter(aufnahmegrund == gr) |> pull(laktat_mmol_l) |> na.omit()
  ci <- t.test(g, conf.level = 0.95)$conf.int
  cat(sprintf("%s: n=%d, Mittel=%.2f, 95%%-KI=[%.2f, %.2f]\n",
              gr, length(g), mean(g), ci[1], ci[2]))
}

from statsmodels.stats.multitest import multipletests

markers = ["leukozyten_g_l", "haemoglobin_g_dl", "kreatinin_mg_dl",
           "laktat_mmol_l", "natrium_mmol_l"]

results = []
for m in markers:
    alive = df.loc[df["verstorben_30d"] == 0, m].dropna()
    dead  = df.loc[df["verstorben_30d"] == 1, m].dropna()
    t_s, p_s = stats.ttest_ind(alive, dead, equal_var=False)
    pooled_sd = np.sqrt(((len(alive)-1)*alive.std()**2 + (len(dead)-1)*dead.std()**2) /
                        (len(alive) + len(dead) - 2))
    d_val = (alive.mean() - dead.mean()) / pooled_sd
    results.append({"marker": m, "t": t_s, "p_raw": p_s, "cohens_d": d_val})

import pandas as pd
res = pd.DataFrame(results)

_, p_bonf, _, _ = multipletests(res["p_raw"], method="bonferroni")
_, p_bh,   _, _ = multipletests(res["p_raw"], method="fdr_bh")
res["p_bonferroni"] = p_bonf
res["p_bh"]         = p_bh
res["sig_bonf"]     = p_bonf < 0.05
res["sig_bh"]       = p_bh   < 0.05

print(res[["marker", "t", "p_raw", "cohens_d", "p_bonferroni", "p_bh",
           "sig_bonf", "sig_bh"]].round(4).to_string(index=False))

suppressPackageStartupMessages(library(tidyverse))
markers <- c("leukozyten_g_l", "haemoglobin_g_dl", "kreatinin_mg_dl",
             "laktat_mmol_l", "natrium_mmol_l")

results <- map_dfr(markers, function(m) {
  alive <- df |> filter(verstorben_30d == 0) |> pull(!!sym(m)) |> na.omit()
  dead  <- df |> filter(verstorben_30d == 1) |> pull(!!sym(m)) |> na.omit()
  test  <- t.test(alive, dead, var.equal = FALSE)
  sp    <- sqrt(((length(alive)-1)*var(alive) + (length(dead)-1)*var(dead)) /
               (length(alive) + length(dead) - 2))
  tibble(marker = m, t = test$statistic, p_raw = test$p.value,
         cohens_d = (mean(alive) - mean(dead)) / sp)
})

results <- results |>
  mutate(p_bonferroni = p.adjust(p_raw, "bonferroni"),
         p_bh         = p.adjust(p_raw, "BH"),
         sig_bonf     = p_bonferroni < 0.05,
         sig_bh       = p_bh < 0.05)
print(results)

from statsmodels.stats.power import TTestIndPower

sepsis    = df.loc[df["aufnahmegrund"] == "Sepsis",  "laktat_mmol_l"].dropna()
no_sepsis = df.loc[df["aufnahmegrund"] != "Sepsis",  "laktat_mmol_l"].dropna()

_, p_t = stats.ttest_ind(sepsis, no_sepsis, equal_var=False)
pooled_sd = np.sqrt(((len(sepsis)-1)*sepsis.std()**2 +
                     (len(no_sepsis)-1)*no_sepsis.std()**2) /
                    (len(sepsis) + len(no_sepsis) - 2))
d_lac = (sepsis.mean() - no_sepsis.mean()) / pooled_sd

analysis = TTestIndPower()
# Observed power given actual sample sizes
power_obs = analysis.power(effect_size=abs(d_lac),
                            nobs1=len(sepsis),
                            ratio=len(no_sepsis) / len(sepsis),
                            alpha=0.05, alternative="two-sided")
# Required n per group for power = 0.8
n_needed = analysis.solve_power(effect_size=abs(d_lac), power=0.8,
                                 ratio=1.0, alpha=0.05)

print(f"Sepsis:    n={len(sepsis)}, mean={sepsis.mean():.2f}")
print(f"kein Sep.: n={len(no_sepsis)}, mean={no_sepsis.mean():.2f}")
print(f"Cohen's d={d_lac:.3f}, t-p={p_t:.4f}")
print(f"Observed power = {power_obs:.3f}")
print(f"n per group for power=0.8: ~{int(np.ceil(n_needed))}")

# install.packages("pwr") # once
library(pwr)

# Parameters from observed data (approximate)
d_lac    <- 0.197   # observed Cohen's d
n_sepsis <- 87
n_other  <- 327

# Harmonic mean n for unequal groups
n_harm <- 2 / (1/n_sepsis + 1/n_other)
power_res <- pwr.t.test(d = d_lac, n = n_harm / 2,
                         sig.level = 0.05, type = "two.sample",
                         alternative = "two.sided")
print(power_res)

# Required n per group for power = 0.8
pwr.t.test(d = d_lac, power = 0.8, sig.level = 0.05,
           type = "two.sample", alternative = "two.sided")