Post-#118 follow-up: merge gompertz scripts + update bias-figure TODO

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.
 
-<!-- TODO: figure needs adaptaiton, I think remove left-truncation, delayed entry and immortal time. these will be dealt with in later chapters. -->
+<!-- TODO (Raphael): restyle figure for visual consistency; keep all 5 rows -->
 
 ![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"}
 

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")