From 65d2943631ea6f16841856a78e4d8f7060b6541a Mon Sep 17 00:00:00 2001 From: adibender Date: Fri, 8 May 2026 15:53:56 +0200 Subject: [PATCH] Merge intro_gompertz_human.R into code.R; update bias-figure TODO --- book/P0C1_intro.qmd | 2 +- book/experiments/code.R | 65 +++++++++++++++++++----- book/experiments/intro_gompertz_human.R | 66 ------------------------- 3 files changed, 54 insertions(+), 79 deletions(-) delete mode 100644 book/experiments/intro_gompertz_human.R diff --git a/book/P0C1_intro.qmd b/book/P0C1_intro.qmd index 2256f39e..53f59b4a 100644 --- a/book/P0C1_intro.qmd +++ b/book/P0C1_intro.qmd @@ -133,7 +133,7 @@ Methods for handling left-truncation, and the detailed treatment of selection an The first part of this book will formally introduce and develop these concepts further. - + ![Common forms of censoring and truncation with biases. First row (fully observed patient): a patient is diagnosed with a disease and immediately enrolled into study and observed to die during the study period. Second row (right censoring): a patient is diagnosed with a disease and immediately enrolled into study but they drop-out before they die, their true survival time is unknown. Third row (left truncation, delayed entry): a patient is diagnosed before they enter the study but enrolled at a later time, they are erroneously treated as if they were diagnosed on the date of enrolment. Fourth row (left truncation selection bias): a patient dies before they enter the study and are never observed. Fifth row (immortal time bias): a patient is diagnosed with a disease and immediately enrolled into study. However, in contrast to the patient in the top row in the control group who can die at any point from enrolment, the patient in the treatment group (bottom row) cannot die between enrolment and treatment (otherwise they could not have entered the treatment group).](Figures/introduction/censoring.png){#fig-intro-censoring fig-alt="TODO"} diff --git a/book/experiments/code.R b/book/experiments/code.R index fe077999..762f14e7 100644 --- a/book/experiments/code.R +++ b/book/experiments/code.R @@ -24,19 +24,60 @@ library(ggpubr) library(pseudo) library(mgcv) library(broom) +library(eha) # swedeaths / swepop + +## Intro - Gompertz fit to real Swedish age-specific mortality (2019, ages 30-100). +## Data source: eha::swedeaths and eha::swepop (Statistics Sweden, bundled in eha). +## Output: book/Figures/introduction/gompertz.png +year_use <- 2019 +age_lo <- 30L +age_hi <- 100L + +d_int <- subset(swedeaths, year == year_use & age >= age_lo & age <= age_hi) +p_int <- subset(swepop, year == year_use & age >= age_lo & age <= age_hi) + +d_agg <- aggregate(deaths ~ age, data = d_int, sum) +p_agg <- aggregate(pop ~ age, data = p_int, sum) +lt <- merge(d_agg, p_agg) +lt$mx <- lt$deaths / lt$pop +lt$x <- lt$age - age_lo + +# Gompertz hazard h(x) = a * exp(b * x), fit by population-weighted log-linear +# regression on the age-specific mortality rate. +fit_gp <- lm(log(mx) ~ x, data = lt, weights = pop) +a_gp <- unname(exp(coef(fit_gp)[1])) +b_gp <- unname(coef(fit_gp)[2]) +cat(sprintf("Gompertz fit: a = %.4g, b = %.4f (age origin = %d)\n", + a_gp, b_gp, age_lo)) + +ages_gp <- seq(age_lo, 110, length.out = 400) +x_gp <- ages_gp - age_lo +hx_gp <- a_gp * exp(b_gp * x_gp) +Hx_gp <- (a_gp / b_gp) * (exp(b_gp * x_gp) - 1) +Sx_gp <- exp(-Hx_gp) +Fx_gp <- 1 - Sx_gp +fx_gp <- hx_gp * Sx_gp + +fun_levels_gp <- c("f(t): density", + "h(t): hazard", + "F(t): CDF", + "S(t): survival") +df_gp <- rbind( + data.frame(age = ages_gp, y = fx_gp, func = "f(t): density"), + data.frame(age = ages_gp, y = hx_gp, func = "h(t): hazard"), + data.frame(age = ages_gp, y = Fx_gp, func = "F(t): CDF"), + data.frame(age = ages_gp, y = Sx_gp, func = "S(t): survival") +) +df_gp$func <- factor(df_gp$func, levels = fun_levels_gp) -## Intro - distributions -g = dstr("Gompertz", shape = 2, decorators = "ExoticStatistics") -t = seq.int(0, 1.5, length.out = 100) -functions <- c("Probability density function, y=f(t)", "Cumulative distribution function, y=F(t)", "Hazard function, y=h(t)", "Survival function, y=S(t)") -d = data.frame(t = t, fun = factor(rep(functions, each = 100), levels = functions), y = c(g$pdf(t), g$hazard(t), g$cdf(t), g$survival(t))) -g = ggplot(d, aes(x = t, y = y, color = fun)) + - geom_line() + - facet_wrap(~fun, scales = "free", nrow = 2) + - theme_bw() + - theme(legend.position = "n") -ggsave("book/Figures/introduction/gompertz.png", g, height = 3, units = "in", - dpi = 600) +g_gp <- ggplot(df_gp, aes(x = age, y = y)) + + geom_line(linewidth = 0.7) + + facet_wrap(~ func, scales = "free_y", nrow = 2) + + theme_bw(base_size = 11) + + labs(x = "age (years)", y = NULL) + +ggsave("book/Figures/introduction/gompertz.png", g_gp, + height = 3.5, width = 7, units = "in", dpi = 600) ## Ranking s_t = tsk("whas") diff --git a/book/experiments/intro_gompertz_human.R b/book/experiments/intro_gompertz_human.R deleted file mode 100644 index a5734bbe..00000000 --- a/book/experiments/intro_gompertz_human.R +++ /dev/null @@ -1,66 +0,0 @@ -# Generate the four-panel Gompertz figure (PDF, hazard, CDF, survival) with -# parameters fit to real Swedish age-specific mortality data from 2019 -# (pre-COVID), both sexes combined, ages 30-100. -# -# Data: eha::swedeaths and eha::swepop -# Source: https://www.scb.se/ (Statistics Sweden), bundled in the R package -# `eha` (covers 1969-2020). -# -# Output: book/Figures/introduction/gompertz.png - -suppressPackageStartupMessages({ - library(eha) - library(ggplot2) -}) - -year_use <- 2019 -age_lo <- 30L -age_hi <- 100L - -d <- subset(swedeaths, year == year_use & age >= age_lo & age <= age_hi) -p <- subset(swepop, year == year_use & age >= age_lo & age <= age_hi) - -d_agg <- aggregate(deaths ~ age, data = d, sum) -p_agg <- aggregate(pop ~ age, data = p, sum) -lt <- merge(d_agg, p_agg) -lt$mx <- lt$deaths / lt$pop -lt$x <- lt$age - age_lo - -## Fit Gompertz hazard h(x) = a * exp(b * x) by population-weighted -## log-linear regression on the age-specific mortality rate. -fit <- lm(log(mx) ~ x, data = lt, weights = pop) -a <- unname(exp(coef(fit)[1])) -b <- unname(coef(fit)[2]) - -cat(sprintf("Gompertz fit: a = %.4g, b = %.4f (age origin = %d)\n", - a, b, age_lo)) - -ages <- seq(age_lo, 110, length.out = 400) -x <- ages - age_lo -hx <- a * exp(b * x) -Hx <- (a / b) * (exp(b * x) - 1) -Sx <- exp(-Hx) -Fx <- 1 - Sx -fx <- hx * Sx - -fun_levels <- c("f(t): density", - "h(t): hazard", - "F(t): CDF", - "S(t): survival") - -df <- rbind( - data.frame(age = ages, y = fx, func = "f(t): density"), - data.frame(age = ages, y = hx, func = "h(t): hazard"), - data.frame(age = ages, y = Fx, func = "F(t): CDF"), - data.frame(age = ages, y = Sx, func = "S(t): survival") -) -df$func <- factor(df$func, levels = fun_levels) - -g <- ggplot(df, aes(x = age, y = y)) + - geom_line(linewidth = 0.7) + - facet_wrap(~ func, scales = "free_y", nrow = 2) + - theme_bw(base_size = 11) + - labs(x = "age (years)", y = NULL) - -ggsave("book/Figures/introduction/gompertz.png", g, - height = 3.5, width = 7, units = "in", dpi = 600)