Motivation

Deaths in older people from falls on the increase ONS figures Independent article

Data:

• Deaths in 21st century https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/datasets/the21stcenturymortalityfilesdeathsdataset/current/regdeaths2001to2016.xls

• GBD http://ihmeuw.org/4e0y

GBD suggests death from falls increased from 2600 in over 70s in 1990 to 4873 in 2016, and death rate from 50 per 100,000 to 73.2.

Note GBD figures are modelled estimates (more later). The charts (see below) suggest increases happened since 2000

library(downloader)
library(readxl)
library(tidyverse, quietly = TRUE)

## ── Attaching packages ──────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

## ✔ ggplot2 2.2.1     ✔ purrr   0.2.4
## ✔ tibble  1.4.2     ✔ dplyr   0.7.4
## ✔ tidyr   0.8.0     ✔ stringr 1.3.0
## ✔ readr   1.1.1     ✔ forcats 0.3.0

## ── Conflicts ─────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(dsrTest)

## Download deaths data
falls <- download("https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/datasets/the21stcenturymortalityfilesdeathsdataset/current/regdeaths2001to2016.xls", "falls.xls")

## Download population data
pop <- download("https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/the21stcenturymortalityfilespopulationdataset/current/populations20012016.xls", "pop.xls")

## Identify sheets
sheets <- readxl::excel_sheets("falls.xls")

## Function to extract sheets to data frames
sheet_extract <- function(file, sheet) {
  
  sheet <- read_excel(file, sheet = sheet, skip = 1)
  

}

## Test
s <- sheet_extract("falls.xls", sheets[3])

## Annoyingly ONS changed headers after 2013 - will import separately
falls_data2001_13 <- map_df(sheets[3:15], function(x) sheet_extract("falls.xls", x))

falls2014 <- sheet_extract("falls.xls", sheets[16])
falls2015_16 <- map_df(sheets[17:18], function(x) sheet_extract("falls.xls", x))

## Convert headers to upper cases and change YEAR to YR
colnames(falls2014) <- toupper(colnames(falls2014))
colnames(falls2015_16) <- toupper(colnames(falls2015_16))
death_14_15 <- bind_rows(falls2014, falls2015_16)
colnames(death_14_15) <- c("ICD-10","YR", "SEX", "AGE", "NDTHS")

## And bind all

deaths_all <- bind_rows(falls_data2001_13, death_14_15)

gbd_falls1 <- read_csv("download (36).csv")

## Parsed with column specification:
## cols(
##   Location = col_character(),
##   Year = col_integer(),
##   Age = col_character(),
##   Sex = col_character(),
##   `Cause of death or injury` = col_character(),
##   Measure = col_character(),
##   Value = col_double(),
##   `Lower bound` = col_double(),
##   `Upper bound` = col_double()
## )

## Warning in rbind(names(probs), probs_f): number of columns of result is not
## a multiple of vector length (arg 1)

## Warning: 3 parsing failures.
## row # A tibble: 3 x 5 col     row col   expected  actual    file                expected   <int> <chr> <chr>     <chr>     <chr>               actual 1    28 <NA>  9 columns 1 columns 'download (36).csv' file 2    29 <NA>  9 columns 1 columns 'download (36).csv' row 3    30 <NA>  9 columns 1 columns 'download (36).csv'

gbd_falls2 <- read_csv("download (37).csv")

## Parsed with column specification:
## cols(
##   Location = col_character(),
##   Year = col_integer(),
##   Age = col_character(),
##   Sex = col_character(),
##   `Cause of death or injury` = col_character(),
##   Measure = col_character(),
##   Value = col_double(),
##   `Lower bound` = col_double(),
##   `Upper bound` = col_double()
## )

## Warning in rbind(names(probs), probs_f): number of columns of result is not
## a multiple of vector length (arg 1)

## Warning: 3 parsing failures.
## row # A tibble: 3 x 5 col     row col   expected  actual    file                expected   <int> <chr> <chr>     <chr>     <chr>               actual 1    28 <NA>  9 columns 1 columns 'download (37).csv' file 2    29 <NA>  9 columns 1 columns 'download (37).csv' row 3    30 <NA>  9 columns 1 columns 'download (37).csv'

gbd_falls <- bind_rows(gbd_falls1, gbd_falls2)

gbd_falls_2001_1670  <- gbd_falls %>%
  filter(Year > 2000) %>%
  mutate(SEX = case_when(Sex == "Male" ~ 1, 
                         Sex == "Female" ~ 2)) %>%
  select(YR = Year, SEX, deaths1 = Value) %>%
  mutate(type = "GBD estimate")

gbd_falls_2001_1670 %>%
  group_by(SEX) %>%
  do(broom::tidy(lm(deaths1~YR, data = .))) %>%
  filter(term == "YR")

Assume ICD10 for falls W00-19 http://www.icd10data.com/ICD10CM/Codes/V00-Y99/W00-W19

falls_deaths <- deaths_all %>%
  filter(str_detect(`ICD-10`, "^W0|^W1"))


x59_deaths70 <- deaths_all %>%
  filter(str_detect(`ICD-10`, "X59")) %>%
  filter(AGE %in% c("70-74", "75-79", "80-84", "85+")) %>%
  group_by(YR, SEX) %>%
  summarise(deaths1 = sum(NDTHS)) %>%
  mutate(type = "x59")

X592001 <- filter(x59_deaths70, YR == 2001, SEX == 2)
X592016 <- filter(x59_deaths70, YR == 2016, SEX == 2)


falls70 <- falls_deaths %>%
  group_by(YR, SEX) %>%
  filter(AGE %in% c("70-74", "75-79", "80-84", "85+")) %>%
  summarise(deaths1 = sum(NDTHS)) %>%
  #mutate(deaths70 = sum(deaths)) %>%
  ungroup() %>%
  mutate(type = "W00-19")

falls_data <- bind_rows(falls70, gbd_falls_2001_1670, x59_deaths70)

falls_data1 <- falls_data %>%
  filter(type != "GBD estimate") %>%
  group_by(YR, SEX) %>%
  summarise(tot = sum(deaths1))

falls_data %>%
  ggplot(aes(YR, deaths1)) +
  geom_line(aes(group = type, colour = type)) +
  geom_line(data = falls_data1, aes(YR, tot), lty = "dotted") +
  facet_wrap(~SEX) +
  geom_vline(xintercept = 2011, lty = "dotted", color = "red") +
  geom_vline(xintercept = 2014, lty = "dotted", color = "green") +
  annotate("text", x = 2011, y = 4000, label= "ICD2001.2 \nto 2010", size = 2.4, family = "Helvetica")+
  annotate("text", x = 2014, y = 4000, label= "IRIS introduced", size = 2.4, family = "Arial")+
  labs(title = "Deaths from falls 2001-2016, 70+", 
       subtitle = "ICD10 W10-19") +
  labs(y = "No deaths") +
  govstyle::theme_gov() +
  theme(legend.position = "bottom") 

garbage <- deaths_all %>%
  mutate(icd10 = substr(`ICD-10`, 1,3)) %>%
  filter(str_detect(`ICD-10`, "^W|^X")) %>%
  arrange(YR, AGE, SEX, icd10) %>%
  filter(!icd10 %in% c("W69", "W70", "W77",
                       "X24", "X26", "X32", 
                       "X33", "X34", "X35",
                       "X36", "X37", "X38", 
                       "X39", "X59", "X91"))

garbage %>%
  group_by(YR, AGE, SEX) %>%
  mutate(tot = sum(NDTHS), 
         cumtot = sum(tot), 
         `%icd` = round(tot/cumtot, 2),
         `%falls` = ifelse(str_detect(icd10, "^W0|^W1"), sum(`%icd`), 0)) %>%
  filter(str_detect(icd10, "^W0|^W1")) %>%
  group_by(YR, AGE, SEX) %>%
  summarise(falls = sum(`%falls`)) %>%
  #spread(SEX, falls, fill = 0) %>%
  filter(AGE %in% c("70-74", "75-79", "80-84", "85+")) %>%
  ggplot(aes(YR, falls, group = AGE, shape = AGE)) +
  geom_point() +
  facet_grid(AGE~SEX) + 
  geom_smooth(se = FALSE, lwd = .4) +
  labs(y = "% garbage codes assigned to falls") +
  govstyle::theme_gov()

## `geom_smooth()` using method = 'loess'