Part 1: Setup.

The first time you run this code, you must “install” some packages. We use the “pacman” package to check for and install all of the other packages we want.

if (!require("pacman")) # If you don't have the required pacman package, this installs it.
  install.packages("pacman") 
## Loading required package: pacman
pacman::p_load(demography, dplyr, ggplot2, tidyr, ggpubr) # Here is our package list. There are over 20,000 R packages to choose from!

Next, please make your own account at https://www.mortality.org/. Click “New User.” A numeric password will be emailed to the email you provide.

Part 2: Loading in data from the Human Mortality Database

# Change the username and password in the command below to your account.
usa <- hmd.mx(country = "USA",
              username = "[enter email here]",
              password = "[enter numeric password here]", 
              label = "United States mortality data")
# You may get an error saying: "NAs introduced by coercion." You can ignore this. 
# Reformat data into a life table, separated by sex, for multiple years. 
out <- list()
for (y in c(1933, 2016)){
  male = lifetable(usa,
                    years=y,
                    ages=usa$age,
                    type=c("period"),
                    series = c("male"))
  female = lifetable(usa,
                   years=y,
                   ages=usa$age,
                   type=c("period"),
                   series = c("female"))
  male_life_table <- bind_cols(male[1:13])
  female_life_table <- bind_cols(female[1:13])
  out[[y]] <- bind_rows(male_life_table, female_life_table) # Combine male and female life tables for year.
}

# Combine data from all years and reformat some columns.
full_lifetable <- bind_rows(out) %>% 
  mutate(year = as.factor(year)) %>%  # Treat year as category instead of continuous
  rename(sex = series) # Rename column to be more informative

# Check the head of the table (top 6 lines) and the dimensions.
head(full_lifetable)
## # A tibble: 6 x 13
##     age year   mx[,1]  qx[,1] lx[,1]  dx[,1] Lx[,1] Tx[,1] ex[,1] rx[,1] sex  
##   <dbl> <fct>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <chr>
## 1     0 1933  0.0682  0.0648   1     0.0648   0.950   59.2   59.2  0.950 male 
## 2     1 1933  0.0100  0.00999  0.935 0.00934  0.931   58.2   62.3  0.980 male 
## 3     2 1933  0.00467 0.00466  0.926 0.00431  0.924   57.3   61.9  0.993 male 
## 4     3 1933  0.00333 0.00333  0.922 0.00307  0.920   56.4   61.2  0.996 male 
## 5     4 1933  0.00254 0.00253  0.919 0.00233  0.917   55.5   60.4  0.997 male 
## 6     5 1933  0.00209 0.00209  0.916 0.00191  0.915   54.6   59.5  0.998 male 
## # … with 2 more variables: type <chr>, label <chr>
dim(full_lifetable)
## [1] 404  13

Part 3: Create figures!

# Plot life expectancy
p1 <- ggplot(full_lifetable) + 
  geom_line(aes(x = age, y = ex, color = year, linetype=sex)) + 
  ylab("ex (life expectancy)") +
  theme_light() +
  theme(plot.caption = element_text(hjust = 0)) 

# Plot survivorship
p2 <- ggplot(full_lifetable) + 
  geom_line(aes(x = age, y = lx, color = year, linetype=sex)) + 
  ylab("lx (survivorship)") +
  theme_light() +
  theme(plot.caption = element_text(hjust = 0)) 

# Plot mortality
p3 <- ggplot(full_lifetable) + 
  geom_line(aes(x = age, y = dx, color = year, linetype=sex)) + 
  ylab("dx (mortality)") +
  theme_light() 

ggarrange(p1, p2, p3, nrow = 1, common.legend = TRUE, legend = "bottom")