In n8thangreen/survivorETHPOP: Survival curves using ETHPOP data

knitr::opts_chunk$set(echo = TRUE)

The raw annual risk of death is calculated by the proportion of people in each age, sex, ethnicity group who die each year (with those who emigrate being censored, of course). Separate survivorship curve for each year of birth, as life expectancy changes over time. There is sex differences; ethnicity doesn't make much of a difference.

ETHPOP is based on ONS data in what year?

library(readr)
library(purrr)
library(dplyr)
library(ggplot2)
library(scales)
library(reshape2)
library(survivorETHPOP)

ETHPOP

We want to compare between the ONS mortality statistics and ETHPOP. From here it details their method. They calculate a central rate of mortality as the average across 3 years. $$ m_x = \sum_{y1,y2,y3} deaths_i/ \sum_{y1,y2,y3} pop_i $$

Finally, they calculate the mortality rate which is what we will be using to compare and is equivalent to the hazard. $$ q_x = 2 m_x/(2 + m_x) $$

Individual categories hazards and survival

ETHPOP_lifetable <- make_ETHPOP_lifetable()
# save(ETHPOP_lifetable, file = here::here("data", "ETHPOP_lifetable.RData"))


head(ETHPOP_lifetable)

ETHPOP_lifetable %>% 
  survivor_curve(group = list(sex = "M",
                              ETH.group = "WHO",
                              year = 2011)) %>% 
  haz_plot()

ETHPOP_lifetable %>% 
  survivor_curve(group = list(sex = "M",
                              ETH.group = "WHO",
                              year = 2049)) %>% 
  haz_plot(add = TRUE)

ETHPOP_lifetable %>% 
  survivor_curve(group = list(sex = "M",
                              ETH.group = "WHO",
                              year = 2011)) %>%
  haz_plot()

for (i in 2012:2043) {
  ETHPOP_lifetable %>% 
    survivor_curve(group = list(sex = "M",
                                ETH.group = "WHO",
                                year = i)) %>% 
    haz_plot(add = TRUE)
}

par(mfrow=c(1,2))
dat <- ETHPOP_lifetable %>% 
  survivor_curve(group = list(sex = "M",
                              ETH.group = "BAN",
                              year = 2011)) %>% 
  haz_plot() %>% 
  surv_plot()

Ethnic groups

ethnic_grps <- c("BAN", "BLA", "BLC", "CHI", "IND", "MIX",
                 "OAS", "OBL", "OTH", "PAK", "WBI", "WHO")

ETHPOP_lifetable %>% 
  survivor_curve(group = list(sex = "M",
                              ETH.group = "BAN",
                              year = 2011)) %>% 
  haz_plot()

for (i in ethnic_grps) {
  ETHPOP_lifetable %>% 
    survivor_curve(group = list(sex = "M",
                                ETH.group = i,
                                year = 2011)) %>% 
    haz_plot(add = TRUE)
}

baseyr_2011 <- ETHPOP_lifetable$id[ETHPOP_lifetable$yr_age == "2011_0"][1]

ETHPOP_lifetable_2011M <- ETHPOP_lifetable %>% 
  filter(id == baseyr_2011,
         sex == "M")

ETHPOP_lifetable_2011F <- ETHPOP_lifetable %>% 
  filter(id == baseyr_2011,
         sex == "F")

## ggplot

ggplot(ETHPOP_lifetable_2011M, aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  # scale_y_continuous(trans='log2') +
  coord_trans(y = "log10") +
  theme_bw()

ggplot(ETHPOP_lifetable_2011M, aes(x = age, y = qx, colour = ETH.group)) +
  geom_line() +
  # scale_y_continuous(trans='log2') +
  coord_trans(y = "log10") +
  theme_bw()

ggplot(ETHPOP_lifetable_2011M, aes(x = age, y = qx, colour = ETH.group)) +
  geom_line() +
  # scale_y_continuous(trans='log2') +
  theme_bw()

ggplot(ETHPOP_lifetable_2011M, aes(x = age, y = S, colour = ETH.group)) +
  geom_line() +
  theme_bw()

---

ONS lifetables

Read in and check.

lifetables <- read_lifetables()
# save(lifetables, file = here::here("data", "lifetables.RData"))


ONS_lifetables <-
  do.call(rbind, lifetables) %>% 
  mutate(new_yr = year < dplyr::lag(year, default = Inf),
         id = cumsum(new_yr)) %>% 
  group_by(id) %>% 
  mutate(baseyr = min(year)) %>% 
  ungroup() %>% 
  select(-id, -new_yr) %>% 
  mutate(baseyr = as.factor(baseyr))

par(mfrow= c(2,3))
for (i in seq_along(lifetables)) {
  plot(x = lifetables[[i]]$year[lifetables[[i]]$sex == "M"],
       y = lifetables[[i]]$qx[lifetables[[i]]$sex == "M"], log = "y", type = "l",
       main = names(lifetables)[i], ylab = "h", xlab = "year")
  lines(x = lifetables[[i]]$year[lifetables[[i]]$sex == "F"],
        y = lifetables[[i]]$qx[lifetables[[i]]$sex == "F"], log = "y", type = "l", col = "red")
}

par(mfrow= c(1,2))
for (j in c("M","F")) {
  plot(x = lifetables[[1]]$age[lifetables[[1]]$sex == j],
       y = lifetables[[1]]$qx[lifetables[[1]]$sex == j], log = "y", type = "l",
       main = j, ylab = "h", xlab = "age")
  for (i in seq_along(lifetables)) {
    lines(x = lifetables[[i]]$age[lifetables[[i]]$sex == j],
          y = lifetables[[i]]$qx[lifetables[[i]]$sex == j], log = "y", type = "l", col = i)
  }
}

ggplot(ONS_lifetables, aes(x = age, y = qx, colour = interaction(sex, baseyr))) +
  geom_line() +
  # scale_y_continuous(trans='log2') +
  coord_trans(y = "log10") +
  theme_bw()

Comparison with ONS and ETHPOP

2011

ggplot(ETHPOP_lifetable_2011M, aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  coord_trans(y = "log10") +
  ggtitle("Male 2011") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2011 & ONS_lifetables$sex == "M", ],
            colour = "black")

ggplot(ETHPOP_lifetable_2011F, aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  ggtitle("Female 2011") +
  coord_trans(y = "log10") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2011 & ONS_lifetables$sex == "F", ],
            colour = "black")

ggplot(ETHPOP_lifetable_2011M, aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  ggtitle("Male 2011") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2011 & ONS_lifetables$sex == "M", ],
            colour = "black") +
  ylim(0, 0.003) + xlim(0, 60)

ggplot(ETHPOP_lifetable_2011F, aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  ggtitle("Female 2011") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2011 & ONS_lifetables$sex == "F", ],
            colour = "black") +
  ylim(0, 0.003) + xlim(0, 60)

2030

baseyr_2030 <- ETHPOP_lifetable$id[ETHPOP_lifetable$yr_age == "2030_0"][1]

ETHPOP_lifetable[ETHPOP_lifetable$id == baseyr_2030 & ETHPOP_lifetable$sex == "M", ] %>% 
  ggplot(aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  coord_trans(y = "log10") +
  ggtitle("Male 2030") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2030 & ONS_lifetables$sex == "M", ],
            colour = "black")

ETHPOP_lifetable[ETHPOP_lifetable$id == baseyr_2030 & ETHPOP_lifetable$sex == "F", ] %>% 
  ggplot(aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  ggtitle("Female 2030") +
  coord_trans(y = "log10") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2030 & ONS_lifetables$sex == "F", ],
            colour = "black")

baseyr_2030 <- ETHPOP_lifetable$id[ETHPOP_lifetable$yr_age == "2030_0"][1]

ETHPOP_lifetable[ETHPOP_lifetable$id == baseyr_2030 & ETHPOP_lifetable$sex == "M", ] %>% 
  ggplot(aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  ggtitle("Male 2030") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2030 & ONS_lifetables$sex == "M", ],
            colour = "black") +
  ylim(0, 0.003) + xlim(0, 40)

ETHPOP_lifetable[ETHPOP_lifetable$id == baseyr_2030 & ETHPOP_lifetable$sex == "F", ] %>% 
  ggplot(aes(x = age, y = death_rate, colour = ETH.group)) +
  geom_line() +
  ggtitle("Female 2030") +
  theme_bw() +
  geom_line(aes(age, qx, colour = "ONS"),
            data = ONS_lifetables[ONS_lifetables$baseyr == 2030 & ONS_lifetables$sex == "F", ],
            colour = "black") +
  ylim(0, 0.003) + xlim(0, 40)

Life expectancy

ONS

ggplot(ONS_lifetables[ONS_lifetables$sex == "M", ], aes(age, ex, colour = baseyr)) +
  geom_line() +
  theme_bw()

ggplot(ONS_lifetables[ONS_lifetables$sex == "F", ], aes(age, ex, colour = baseyr)) +
  geom_line() +
  theme_bw()

rmarkdown::render("vignettes/analysis.Rmd", output_dir = "docs/", output_format = "html_document")

n8thangreen/survivorETHPOP documentation built on May 10, 2020, 6:06 p.m.