#---------------------------------------------------------------------------#
# Nom : Graphs.R.                                                           #
# Description : Code used to generate the graphs                            #
# Author(s) : Pietro Violo and Nadine Ouellette                             #
# Date : June 20 2025.                                                      #
# Modifications :                                                           #
#---------------------------------------------------------------------------#

# Notes:
# - This code generates all figures for the Quebec analysis.
# - To replicate for Ontario, change the input data sources accordingly:
#     - Replace `Quebec_obituaries_extracted` with Ontario's equivalent.
#     - Replace CHMD files with Ontario-specific exposures and deaths data.
# - Each section below corresponds to a figure or table from the paper.
# - Graphs are listed first, followed by tables.

# - Turns off scientific notation
options(scipen=999)

#---------------------------------------------------------------------------#
# Package library                                                           #
#---------------------------------------------------------------------------#

library(tidyverse)
library(stringi)
library(ggpubr)
library(xtable)
library(MortalitySmooth)

#---------------------------------------------------------------------------#
# Data
#---------------------------------------------------------------------------#

###### PLEASE NOTE! ######
# The AI-generated obituary text data used in this analysis is random. 
# As a result, the findings and graphs produced from this data will not match 
# the results reported in the paper.

# --- 1. Clear Workspace ---
rm(list = ls())

# --- 2. Load Quebec Obituaries Data ---
load("./Data/Quebec_obituaries_extracted.RData")

# --- 3. Load CHMD Exposures: Collapse 105+ and reshape ---
chmd_exposures <- read.table("./Data/Exposures_1x1.txt", 
                             skip = 1, header = TRUE) %>%
  mutate(
    Age = ifelse(Age == "110+", "110", Age),
    Age = as.numeric(Age),
    Age = ifelse(Age >= 105, 105, Age)
  ) %>%
  group_by(Year, Age) %>%
  summarise(
    Female = sum(Female),
    Male = sum(Male),
    Total = sum(Total),
    .groups = "drop"
  ) %>%
  filter(Year %in% 2017:2022) %>%
  select(-Total) %>%
  pivot_longer(cols = c(Female, Male), names_to = "Gender", values_to = "Exposures") %>%
  mutate(Gender = recode(Gender, "Male" = "Man", "Female" = "Woman"))

# --- 4. Load CHMD Death Counts: Collapse 105+ and reshape ---
chmd_deaths_raw <- read.csv("./Data/Deaths_1x1.txt", 
                            skip = 2, header = TRUE, sep = "\t")

# If needed, fix malformed import
chmd_deaths <- read.table(textConnection(chmd_deaths_raw[, ]))

colnames(chmd_deaths) <- c("Year", "Age", "Female", "Male", "Total")

chmd_deaths <- chmd_deaths %>%
  mutate(
    Age = ifelse(Age == "110+", "110", Age),
    Age = as.numeric(Age),
    Age = ifelse(Age >= 105, 105, Age)
  ) %>%
  group_by(Year, Age) %>%
  summarise(
    Female = sum(Female),
    Male = sum(Male),
    Total = sum(Total),
    .groups = "drop"
  ) %>%
  filter(Year %in% 2017:2022) %>%
  pivot_longer(cols = c(Female, Male), names_to = "Gender", values_to = "Deaths") %>%
  mutate(Gender = recode(Gender, "Male" = "Man", "Female" = "Woman"))

# --- 5. Preprocess Quebec Obituaries: Collapse 105+ and aggregate ---
obituaries_quebec <- obituaries_quebec %>%
  rename(Age = age_at_death) %>%
  mutate(Age = ifelse(Age >= 105, 105, Age)) %>%
  group_by(Age, Gender, Year) %>%
  summarise(Deaths = n(), Source = "Obituaries", .groups = "drop") %>%
  group_by(Gender, Year) %>%
  mutate(Total = sum(Deaths))



###############################################################################
# Figure 1 : Observed death counts per age
###############################################################################

# --- 1. Prepare CHMD Data: Calculate % of deaths per group ---
chmd_pct <- chmd_deaths %>%
  group_by(Year, Gender) %>%
  mutate(
    Total = sum(Deaths),
    pct = Deaths / Total * 100
  ) %>%
  ungroup() %>%
  mutate(Source = "CHMD") %>%
  select(Year, Age, Gender, Deaths, Total, pct, Source)

# --- 2. Combine CHMD and Obituary Data ---
df <- bind_rows(
  obituaries_quebec %>% select(Year, Age, Gender, Deaths, Total, Source),
  chmd_pct %>% select(Year, Age, Gender, Deaths, Total, Source)
) %>%
  mutate(
    Year = as.numeric(Year),
    Age = as.numeric(Age),
    Deaths = as.numeric(Deaths),
    Total = as.numeric(Total)
  )

# --- 3. Aggregate Death Counts for 2017–2022 ---
df20172022 <- df %>%
  filter(Year %in% 2017:2022) %>%
  group_by(Age, Gender, Source) %>%
  summarise(Deaths = sum(Deaths), .groups = "drop") %>%
  group_by(Gender, Source) %>%
  mutate(
    Total = sum(Deaths),
    pct = Deaths / Total * 100
  )

# --- 4. Create Plot: Death pyramid (absolute death counts by age and gender) ---
figure1 <- df20172022 %>%
  filter(Gender == "Woman", Source == "CHMD") %>%
  ggplot(aes(x = Age, y = Deaths)) +
  
  geom_bar(stat = "identity", fill = "#9d0208") +
  
  geom_bar(
    data = df20172022 %>% filter(Gender == "Man", Source == "CHMD"),
    aes(x = Age, y = -Deaths),
    stat = "identity", fill = "navy"
  ) +
  
  geom_bar(
    data = df20172022 %>% filter(Gender == "Woman", Source == "Obituaries"),
    aes(x = Age, y = Deaths),
    stat = "identity", fill = "#ff758f", inherit.aes = FALSE
  ) +
  
  geom_bar(
    data = df20172022 %>% filter(Gender == "Man", Source == "Obituaries"),
    aes(x = Age, y = -Deaths),
    stat = "identity", fill = "cornflowerblue", alpha = 0.7
  ) +
  
  labs(
    y = "Death count",
    x = "Age",
    title = ""
  ) +
  
  coord_flip() +
  
  scale_x_continuous(
    breaks = seq(0, 105, 5),
    labels = c(seq(0, 100, 5), "105+")
  ) +
  
  scale_y_continuous(
    breaks = seq(-8750, 8750, 1750),
    labels = abs(seq(-8750, 8750, 1750)),
    limits = c(-8750, 8750)
  ) +
  
  ggthemes::theme_fivethirtyeight() +
  theme(
    axis.title.x = element_text(),
    axis.title.y = element_text(),
    panel.background = element_rect(fill = "white", colour = "white"),
    plot.background = element_rect(fill = "white", colour = "white")
  )


# Save file
ggsave(plot = figure1,
      filename = "./Outputs/figure1.png",
      height = 2600, width =2700, units = "px", dpi = 600,
      device = "png",limitsize = FALSE,bg="white")


###############################################################################
# Figure 2 : Obituary coverage
###############################################################################

# --- 1. Prepare Aggregated Coverage Data (overall, not plotted here) ---
total_coverage <- df20172022 %>%
  select(Age, Gender, Source, Deaths) %>%
  group_by(Age, Source) %>%
  summarise(Deaths = sum(Deaths), .groups = "drop") %>%
  pivot_wider(names_from = Source, values_from = Deaths) %>%
  mutate(Coverage = Obituaries / CHMD * 100)

# --- 2. Create Plot: Obituary coverage (%) by age and gender ---
figure2 <- df20172022 %>%
  select(Age, Gender, Source, Deaths) %>%
  pivot_wider(names_from = Source, values_from = Deaths) %>%
  mutate(Coverage = Obituaries / CHMD * 100) %>%
  ggplot(aes(x = Age, y = Coverage, color = Gender, group = Gender)) +
  geom_line(size = 1) +
  labs(
    y = "Obituary coverage (%)",
    title = ""
  ) +
  scale_y_continuous(
    breaks = seq(0, 100, 25),
    labels = seq(0, 100, 25),
    limits = c(0, 100)
  ) +
  scale_x_continuous(
    breaks = seq(0, 105, 5),
    labels = c(seq(0, 100, 5), "105+")
  ) +
  scale_color_manual(values = c("Man" = "cornflowerblue", "Woman" = "#ff758f")) +
  ggthemes::theme_fivethirtyeight() +
  theme(
    axis.title.x = element_text(),
    axis.title.y = element_text(),
    panel.background = element_rect(fill = "white", colour = "white"),
    plot.background = element_rect(fill = "white", colour = "white"),
    legend.position = "none"
  )


# Save plot
ggsave(plot = figure2,
      filename = "./Outputs/figure2.png",
      height = 2000, width =3600, units = "px", dpi = 600,
      device = "png",limitsize = FALSE,bg="white")



###############################################################################
# Figure A-1 and A-2 : Obituary coverage per year
###############################################################################

# --- 1. Prepare Yearly Deaths Data by Gender and Source ---
df20172022_year <- df %>%
  filter(Year %in% 2017:2022) %>%
  group_by(Age, Gender, Source, Year) %>%
  summarise(Deaths = sum(Deaths), .groups = "drop") %>%
  group_by(Gender, Source) %>%
  mutate(
    Total = sum(Deaths, na.rm = TRUE),
    pct = Deaths / Total * 100
  )

# --- 2. Compute Overall Coverage per Year and Age ---
total_coverage_year <- df20172022_year %>%
  select(Age, Gender, Source, Deaths, Year) %>%
  group_by(Age, Source, Year) %>%
  summarise(Deaths = sum(Deaths), .groups = "drop") %>%
  pivot_wider(names_from = Source, values_from = Deaths) %>%
  mutate(
    Obituaries = ifelse(is.na(Obituaries), 0, Obituaries),
    Coverage = Obituaries / CHMD * 100
  )

# --- 3. Define Plotting Function for a Given Year ---
make_coverage_plot <- function(df, year) {
  df %>%
    filter(Year == year) %>%
    select(Age, Gender, Source, Deaths, Year) %>%
    pivot_wider(names_from = Source, values_from = Deaths) %>%
    mutate(
      Obituaries = ifelse(is.na(Obituaries), 0, Obituaries),
      Coverage = Obituaries / CHMD * 100
    ) %>%
    ggplot(aes(x = Age, y = Coverage, color = Gender, group = Gender)) +
    geom_line(size = 1) +
    scale_y_continuous(
      breaks = seq(0, 100, 25),
      labels = seq(0, 100, 25),
      limits = c(0, 100)
    ) +
    scale_x_continuous(
      breaks = seq(0, 105, 5),
      labels = c(seq(0, 100, 5), "105+")
    ) +
    scale_color_manual(values = c("Man" = "cornflowerblue", "Woman" = "#ff758f")) +
    labs(
      y = "Obituary coverage (%)",
      title = as.character(year)
    ) +
    ggthemes::theme_fivethirtyeight() +
    theme(
      panel.background = element_rect(fill = "white", colour = "white"),
      plot.background  = element_rect(fill = "white", colour = "white"),
      axis.title.x = element_text(),
      axis.title.y = element_text(),
      legend.position = "none"
    )
}

# --- 4. Generate Yearly Plots ---
years <- 2017:2022
plots <- purrr::map(years, ~ make_coverage_plot(df20172022_year, .x))

# --- 5. Combine Plots into One Figure ---
figure_a1a2 <- ggpubr::ggarrange(plotlist = plots, nrow = 6)

# Save Figure
ggsave(plot = figure_a1a2,
      filename = "./Outputs/figurea1_a2.png",
      height = 7200, width = 5280, units = "px", dpi = 600,
      device = "png", limitsize = FALSE, bg = "white")



###############################################################################
# Figure 3 : Proportion of women's deaths
###############################################################################

# --- 1. Prepare Data: Compute proportion of deaths that are women ---
gender_dist_totals <- df20172022 %>%
  select(-Total, -pct) %>%
  pivot_wider(names_from = Gender, values_from = Deaths) %>%
  mutate(
    Total = Man + Woman,
    Woman = Woman / Total
  )

# --- 2. Create Plot: Proportion of women's deaths by age and source ---
figure3 <- gender_dist_totals %>%
  ggplot(aes(x = Age, y = Woman * 100, color = Source, group = Source)) +
  geom_line(size = 1) +
  lims(x = c(0, 105), y = c(0, 100)) +
  scale_x_continuous(
    breaks = seq(0, 105, 5),
    labels = c(seq(0, 100, 5), "105+")
  ) +
  scale_color_manual(values = c("CHMD" = "#9d0208", "Obituaries" = "hotpink")) +
  ggthemes::theme_fivethirtyeight() +
  theme(
    axis.title.x = element_text(),
    axis.title.y = element_text(),
    panel.background = element_rect(fill = "white", colour = "white"),
    plot.background = element_rect(fill = "white", colour = "white")
  ) +
  labs(
    y = "Proportion (%)",
    x = "Age"
  )


# Save plot
ggsave(plot = figure3,
      filename = "./Outputs/figure3.png",
      height = 2000, width =3600, units = "px", dpi = 600,
      device = "png",limitsize = FALSE,bg="white")


###############################################################################
# Figure A3-A4 : Proportion of women's deaths
###############################################################################

# --- 1. Prepare data: Calculate proportion of deaths that are women by year ---
df20172022_year <- df %>%
  filter(Year %in% 2017:2022) %>%
  group_by(Year, Age, Gender, Source) %>%
  summarise(Deaths = sum(Deaths), .groups = "drop") %>%
  pivot_wider(
    names_from = Gender,
    values_from = Deaths,
    values_fill = list(Man = 0, Woman = 0)
  ) %>%
  mutate(
    pct_women = Woman / (Woman + Man) * 100
  )

# --- 2. Define function to plot % of deaths that are women by age ---
make_coverage_plot <- function(data, year) {
  data_year <- data %>% filter(Year == year)
  
  ggplot(data_year, aes(x = Age + 0.5, y = pct_women, color = Source)) +
    geom_line(linewidth = 1.3) +
    labs(
      title = paste0(year),
      x = "Age",
      y = "% of deaths that are women"
    ) +
    scale_color_manual(
      values = c(
        "CHMD" = "#9d0208",
        "Obituaries" = "hotpink"
      ),
      labels = c(
        "CHMD" = "Woman.CHMD",
        "Obituaries" = "Woman.Obits"
      ),
      name = NULL
    ) +
    ggthemes::theme_fivethirtyeight() +
    theme(
      axis.title.y = element_text(),
      axis.title.x = element_text(),
      panel.background = element_rect(fill = 'white', colour = 'white'),
      plot.background = element_rect(fill = 'white', colour = 'white'),
      legend.position = "bottom"
    ) +
    ylim(0, 100) +
    scale_x_continuous(
      breaks = seq(0, 105, 5),
      labels = c(seq(0, 100, 5), "105+")
    )
}

# --- 3. Generate one plot per year ---
years <- 2017:2022
plots <- map(years, ~ make_coverage_plot(df20172022_year, .x))

# --- 4. Combine plots into one figure ---
figurea3a4 <- do.call(ggarrange, c(plots, nrow = 6))


# Save file
ggsave(
  plot = figurea3a4,
  filename = "./Outputs/figurea3a4.png",
  height = 6000, width = 3900, units = "px", dpi = 500,
  device = "png", limitsize = FALSE, bg = "white"
)



###############################################################################
# Figure 4: Age Distribution of Observed Deaths (2017–2022)
###############################################################################

# --- 1. Create Plot: Proportion of deaths by age, gender, and source ---
figure4 <- df20172022 %>%
  ggplot(aes(x = Age + 0.5, y = pct)) +
  
  # Male – Obituaries
  geom_line(
    data = filter(df20172022, Gender == "Man", Source == "Obituaries"),
    linewidth = 1.3,
    color = "cornflowerblue"
  ) +
  
  # Male – CHMD
  geom_line(
    data = filter(df20172022, Gender == "Man", Source == "CHMD"),
    aes(x = Age + 0.5, y = pct),
    linewidth = 1.3,
    color = "navy"
  ) +
  
  # Female – Obituaries
  geom_line(
    data = filter(df20172022, Gender == "Woman", Source == "Obituaries"),
    aes(x = Age + 0.5, y = pct),
    linewidth = 1.3,
    color = "hotpink"
  ) +
  
  # Female – CHMD
  geom_line(
    data = filter(df20172022, Gender == "Woman", Source == "CHMD"),
    aes(x = Age + 0.5, y = pct),
    linewidth = 1.3,
    color = "#9d0208"
  ) +
  
  labs(
    title = "",
    y = "Proportion of deaths (%)",
    x = "Age"
  ) +
  
  scale_x_continuous(
    breaks = seq(0, 105, 5),
    labels = c(seq(0, 100, 5), "105+")
  ) +
  
  coord_cartesian(ylim = c(0, 5)) +
  
  ggthemes::theme_fivethirtyeight() +
  theme(
    axis.title.y = element_text(),
    axis.title.x = element_text(),
    panel.background = element_rect(fill = "white"),
    plot.background  = element_rect(fill = "white")
  )

# Save Figure
ggsave(
  plot = figure4,
  filename = "./Outputs/figure4.png",
  height = 2000, width = 3600, units = "px", dpi = 600,
  device = "png", bg = "white", limitsize = FALSE
)


###############################################################################
# Figures A5–A6: Annual Age Distribution of Observed Deaths (2017–2022)
###############################################################################

# --- 1. Prepare data: Calculate proportion of deaths by year, gender, and source ---
df20172022_year <- df %>%
  filter(Year %in% 2017:2022) %>%
  group_by(Age, Gender, Source, Year) %>%
  summarise(Deaths = sum(Deaths, na.rm = TRUE), .groups = "drop") %>%
  group_by(Gender, Source, Year) %>%
  mutate(
    Total = sum(Deaths, na.rm = TRUE),
    pct = Deaths / Total * 100
  ) %>%
  filter(!is.na(Gender), !is.na(Source)) %>%
  ungroup()

# --- 2. Define function to plot proportion of deaths by age ---
create_year_plot <- function(data, year) {
  df_year <- filter(data, Year == year)
  
  ggplot(df_year, aes(x = Age + 0.5, y = pct, color = interaction(Gender, Source))) +
    geom_line(linewidth = 1.3) +
    labs(
      title = paste(year),
      y = "Proportion of deaths (%)",
      x = "Age"
    ) +
    scale_color_manual(
      values = c(
        "Man.CHMD"          = "navy",
        "Man.Obituaries"    = "cornflowerblue",
        "Woman.CHMD"        = "#9d0208",
        "Woman.Obituaries"  = "hotpink"
      ),
      labels = c(
        "Man.CHMD"          = "Men – CHMD",
        "Man.Obituaries"    = "Men – Obituaries",
        "Woman.CHMD"        = "Women – CHMD",
        "Woman.Obituaries"  = "Women – Obituaries"
      ),
      name = NULL
    ) +
    scale_x_continuous(
      breaks = seq(0, 105, 5),
      labels = c(seq(0, 100, 5), "105+")
    ) +
    coord_cartesian(ylim = c(0, 5)) +
    theme(
      axis.title.x = element_text(),
      axis.title.y = element_text(),
      legend.position = "bottom"
    )
}


# --- 4. Combine plots into one figure ---
figurea5a6 <- ggpubr::ggarrange(create_year_plot(df20172022_year, 2017),
                                    create_year_plot(df20172022_year, 2018),
                                    create_year_plot(df20172022_year, 2019),
                                    create_year_plot(df20172022_year, 2020),
                                    create_year_plot(df20172022_year, 2021),
                                    create_year_plot(df20172022_year, 2022),
                                    nrow = 6)

# --- 5. Save figure to file ---
ggsave(
  plot = figurea5a6,
  filename = "./Outputs/figure_a5a6.png",
  height = 6000, width = 3900, units = "px", dpi = 600,
  device = "png", bg = "white", limitsize = FALSE
)


###############################################################################
# Figure 5-6
###############################################################################

# --- 1. Merge CHMD and Obituary Data: Adjust obituary deaths to CHMD totals ---
df_combined <- left_join(
  chmd_deaths %>% group_by(Gender, Year) %>% mutate(Total = sum(Deaths)),
  obituaries_quebec %>% group_by(Gender, Year) %>%
    mutate(Total = sum(Deaths), pct = Deaths / Total) %>%
    ungroup() %>%
    select(Age, Gender, Year, pct)
) %>%
  mutate(pct = ifelse(is.na(pct), 0, pct),
         Deaths_obit = Total * pct)

# --- 2. Attach Exposure Data to Combined Dataset ---
df_combined <- left_join(df_combined, chmd_exposures)

# --- 3. Compute Life Tables for CHMD (HMD Protocol V6) ---
chmd_life_tables <- df_combined %>%
  select(-Deaths_obit) %>%
  mutate(
    m_x = Deaths / Exposures,  # central death rate
    a_x = case_when(
      Age == 0 & Gender == "Man"   & (m_x <= 0.02300) ~ 0.14929 - 1.99545 * m_x,
      Age == 0 & Gender == "Man"   & (m_x <= 0.08307) ~ 0.02832 + 3.26021 * m_x,
      Age == 0 & Gender == "Man"   & (m_x >= 0.08307) ~ 0.29915,
      Age == 0 & Gender == "Woman" & (m_x <= 0.01724) ~ 0.14903 - 2.05527 * m_x,
      Age == 0 & Gender == "Woman" & (m_x <= 0.06891) ~ 0.04667 + 3.88089 * m_x,
      Age == 0 & Gender == "Woman" & (m_x >= 0.06891) ~ 0.31411,
      Age == 105 ~ ifelse(m_x == 0, 0, 1 / m_x),
      TRUE ~ 1 / 2
    ),
    q_x = case_when(Age == 105 ~ 1, TRUE ~ m_x / (1 + (1 - a_x) * m_x)),
    p_x = 1 - q_x
  ) %>%
  group_by(Gender, Year) %>%
  arrange(Gender, Year) %>%
  mutate(
    l_x = case_when(Age == 0 ~ 100000, TRUE ~ 100000 * lag(cumprod(p_x))),
    d_x = ifelse(Age == 105, l_x, l_x * q_x),
    L_x = ifelse(Age == 105, l_x * a_x, l_x - (1 - a_x) * d_x)
  ) %>%
  arrange(Gender, Year, desc(Age)) %>%
  mutate(T_x = cumsum(L_x), e_x = T_x / l_x) %>%
  arrange(Gender, Year, Age) %>%
  ungroup()

# --- 4. Compute Life Tables for Obituaries ---
obits_life_tables <- df_combined %>%
  select(-Deaths) %>%
  mutate(
    Deaths_obit = ifelse(is.na(Deaths_obit), 0, Deaths_obit),
    m_x = Deaths_obit / Exposures,
    a_x = case_when(
      Age == 0 & Gender == "Man"   & (m_x <= 0.02300) ~ 0.14929 - 1.99545 * m_x,
      Age == 0 & Gender == "Man"   & (m_x <= 0.08307) ~ 0.02832 + 3.26021 * m_x,
      Age == 0 & Gender == "Man"   & (m_x >= 0.08307) ~ 0.29915,
      Age == 0 & Gender == "Woman" & (m_x <= 0.01724) ~ 0.14903 - 2.05527 * m_x,
      Age == 0 & Gender == "Woman" & (m_x <= 0.06891) ~ 0.04667 + 3.88089 * m_x,
      Age == 0 & Gender == "Woman" & (m_x >= 0.06891) ~ 0.31411,
      Age == 105 ~ ifelse(m_x == 0 | is.na(m_x), 0, 1 / m_x),
      TRUE ~ 1 / 2
    ),
    q_x = case_when(Age == 105 ~ 1, TRUE ~ m_x / (1 + (1 - a_x) * m_x)),
    p_x = 1 - q_x
  ) %>%
  group_by(Gender, Year) %>%
  arrange(Gender, Year) %>%
  mutate(
    l_x = case_when(Age == 0 ~ 100000, TRUE ~ 100000 * lag(cumprod(p_x))),
    d_x = ifelse(Age == 105, l_x, l_x * q_x),
    L_x = ifelse(Age == 105, l_x * a_x, l_x - (1 - a_x) * d_x)
  ) %>%
  arrange(Gender, Year, desc(Age)) %>%
  mutate(T_x = cumsum(L_x), e_x = T_x / l_x) %>%
  arrange(Gender, Year, Age) %>%
  ungroup()

# --- 5. Aggregate 2017–2022 Data for Combined Life Tables ---
df_combined_20172022 <- df_combined %>%
  group_by(Age, Gender) %>%
  summarise(
    Deaths = sum(Deaths),
    Deaths_obit = sum(Deaths_obit),
    Exposures = sum(Exposures),
    .groups = "drop"
  )

# --- 6. CHMD Life Table (2017–2022 Combined) ---
chmd_life_tables_20172022 <- df_combined_20172022 %>%
  select(-Deaths_obit) %>%
  mutate(
    m_x = Deaths / Exposures,
    a_x = case_when(
      Age == 0 & Gender == "Man"   & (m_x <= 0.02300) ~ 0.14929 - 1.99545 * m_x,
      Age == 0 & Gender == "Man"   & (m_x <= 0.08307) ~ 0.02832 + 3.26021 * m_x,
      Age == 0 & Gender == "Man"   & (m_x >= 0.08307) ~ 0.29915,
      Age == 0 & Gender == "Woman" & (m_x <= 0.01724) ~ 0.14903 - 2.05527 * m_x,
      Age == 0 & Gender == "Woman" & (m_x <= 0.06891) ~ 0.04667 + 3.88089 * m_x,
      Age == 0 & Gender == "Woman" & (m_x >= 0.06891) ~ 0.31411,
      Age == 105 ~ ifelse(m_x == 0 | is.na(m_x), 0, 1 / m_x),
      TRUE ~ 1 / 2
    ),
    q_x = case_when(Age == 105 ~ 1, TRUE ~ m_x / (1 + (1 - a_x) * m_x)),
    p_x = 1 - q_x
  ) %>%
  group_by(Gender) %>%
  arrange(Gender) %>%
  mutate(
    l_x = case_when(Age == 0 ~ 100000, TRUE ~ 100000 * lag(cumprod(p_x))),
    d_x = ifelse(Age == 105, l_x, l_x * q_x),
    L_x = ifelse(Age == 105, l_x * a_x, l_x - (1 - a_x) * d_x)
  ) %>%
  arrange(Gender, desc(Age)) %>%
  mutate(T_x = cumsum(L_x), e_x = T_x / l_x) %>%
  arrange(Gender, Age) %>%
  ungroup()

# --- 7. Obituary Life Table (2017–2022 Combined) ---
obits_life_tables_20172022 <- df_combined_20172022 %>%
  select(-Deaths) %>%
  mutate(
    m_x = Deaths_obit / Exposures,
    a_x = case_when(
      Age == 0 & Gender == "Man"   & (m_x <= 0.02300) ~ 0.14929 - 1.99545 * m_x,
      Age == 0 & Gender == "Man"   & (m_x <= 0.08307) ~ 0.02832 + 3.26021 * m_x,
      Age == 0 & Gender == "Man"   & (m_x >= 0.08307) ~ 0.29915,
      Age == 0 & Gender == "Woman" & (m_x <= 0.01724) ~ 0.14903 - 2.05527 * m_x,
      Age == 0 & Gender == "Woman" & (m_x <= 0.06891) ~ 0.04667 + 3.88089 * m_x,
      Age == 0 & Gender == "Woman" & (m_x >= 0.06891) ~ 0.31411,
      Age == 105 ~ ifelse(m_x == 0 | is.na(m_x), 0, 1 / m_x),
      TRUE ~ 1 / 2
    ),
    q_x = case_when(Age == 105 ~ 1, TRUE ~ m_x / (1 + (1 - a_x) * m_x)),
    p_x = 1 - q_x
  ) %>%
  group_by(Gender) %>%
  arrange(Gender) %>%
  mutate(
    l_x = case_when(Age == 0 ~ 100000, TRUE ~ 100000 * lag(cumprod(p_x))),
    d_x = ifelse(Age == 105, l_x, l_x * q_x),
    L_x = ifelse(Age == 105, l_x * a_x, l_x - (1 - a_x) * d_x)
  ) %>%
  arrange(Gender, desc(Age)) %>%
  mutate(T_x = cumsum(L_x), e_x = T_x / l_x) %>%
  arrange(Gender, Age) %>%
  ungroup()

# --- 8. Merge Life Table Results for Graphing ---
df_20172022 <- bind_rows(
  obits_life_tables_20172022 %>% mutate(Source = "Obituaries") %>% rename(Deaths = Deaths_obit),
  chmd_life_tables_20172022 %>% mutate(Source = "CHMD")
) %>%
  mutate(Source = paste(Gender, Source))

# --- 9. Plot A: Mortality Rate (q_x, log scale) ---
qx_combined <- df_20172022 %>%
  ggplot(aes(x = Age + 0.5, y = q_x, color = Source, group = Source)) +
  labs(
    y = "Death rate (log scale)",
    color = "Source of data",
    title = "A. Mortality rate"
  ) +
  ggthemes::theme_fivethirtyeight() +
  geom_line(linewidth = 0.8) +
  scale_color_manual(values = c("navy", "cornflowerblue", "#9d0208", "hotpink")) +
  scale_y_log10(limits = c(1e-5, 1)) +
  scale_x_continuous(breaks = seq(0, 105, 5), labels = c(seq(0, 100, 5), "105+")) +
  theme(
    panel.background = element_rect(fill = 'white', colour = 'white'),
    plot.background  = element_rect(fill = 'white', colour = 'white'),
    legend.background = element_rect(fill = 'white', colour = 'white'),
    legend.key = element_rect(fill = 'white', colour = 'white'),
    axis.title.y = element_text()
  ) +
  guides(color = guide_legend(nrow = 2))

# --- 10. Plot B: Survival Function (l_x) ---
lx_combined <- df_20172022 %>%
  ggplot(aes(x = Age + 0.5, y = l_x / 100000, color = Source, group = Source)) +
  labs(
    y = "Probability of survival",
    color = "Source of data",
    title = "B. Survival function"
  ) +
  ggthemes::theme_fivethirtyeight() +
  geom_line(linewidth = 0.8) +
  scale_color_manual(values = c("navy", "cornflowerblue", "#9d0208", "hotpink")) +
  scale_x_continuous(breaks = seq(0, 105, 5), labels = c(seq(0, 100, 5), "105+")) +
  scale_y_continuous(breaks = seq(0, 1, 0.2), labels = seq(0, 1, 0.2), limits = c(0, 1)) +
  theme(
    panel.background = element_rect(fill = 'white', colour = 'white'),
    plot.background  = element_rect(fill = 'white', colour = 'white'),
    legend.background = element_rect(fill = 'white', colour = 'white'),
    legend.key = element_rect(fill = 'white', colour = 'white'),
    axis.title.y = element_text()
  ) +
  guides(color = guide_legend(nrow = 2))

# --- 11. Plot C: Distribution of Deaths (d_x) ---
dx_combined <- df_20172022 %>%
  ggplot(aes(x = Age + 0.5, y = d_x, color = Source, group = Source)) +
  labs(
    y = "Number of life table deaths",
    x = "Age",
    color = "Source of data",
    title = "C. Distribution of deaths"
  ) +
  ggthemes::theme_fivethirtyeight() +
  geom_line(linewidth = 0.8) +
  scale_color_manual(values = c("navy", "cornflowerblue", "#9d0208", "hotpink")) +
  scale_x_continuous(breaks = seq(0, 105, 5), labels = c(seq(0, 100, 5), "105+")) +
  theme(
    panel.background = element_rect(fill = 'white', colour = 'white'),
    plot.background  = element_rect(fill = 'white', colour = 'white'),
    legend.background = element_rect(fill = 'white', colour = 'white'),
    legend.key = element_rect(fill = 'white', colour = 'white'),
    axis.title = element_text()
  ) +
  guides(color = guide_legend(nrow = 2))

# --- 12. Combine and Export Figure Panels ---
combined_plot <- ggpubr::ggarrange(
  qx_combined,
  lx_combined,
  dx_combined,
  common.legend = TRUE,
  nrow = 3,
  legend = "bottom",
  align = "v"
)

# Save figure
ggsave(
  plot = combined_plot,
  filename = "./Outputs/figure5-6.png",
  height = 5000, width = 3900, units = "px", dpi = 600,
  device = "png", limitsize = FALSE, bg = "white"
)



###############################################################################
# Tables
###############################################################################

# Tables 1 and 2 were generated by randomly sampling 1,600 observations 
# distributed across four age groups, two genders, and two provinces. 
# We manually labeled age and gender, then compared these labels with the 
# outputs produced by our automated extraction methods. 
# Since the sample is randomized each time and the obituary texts 
# cannot be shared due to privacy concerns, the replication code 
# begins with Table 3 and onward.

###############################################################################
# Table 3
###############################################################################

# --- Reset environment, keeping only required tables ---
rm(list = setdiff(ls(), c("chmd_deaths", "chmd_exposures", 
                          "chmd_life_tables", "obits_life_tables")))

# --- 1. Load obituary data ---
load("./Data/Quebec_obituaries_extracted.RData")
load("./Data/Ontario_obituaries_extracted.RData")

# --- 2. Helper function to count subsets ---
count_subsets <- function(df, lang) {
  list(
    nrow(df %>% filter(language == lang)),
    nrow(df %>% filter(language == lang, !is.na(Gender) & !is.na(age_at_death))),
    nrow(df %>% filter(language == lang, Gender == "Woman", !is.na(Gender) & !is.na(age_at_death))),
    nrow(df %>% filter(language == lang, Gender == "Man",   !is.na(Gender) & !is.na(age_at_death))),
    nrow(df %>% filter(language == lang, !is.na(Gender) &  is.na(age_at_death))),
    nrow(df %>% filter(language == lang, Gender == "Woman", !is.na(Gender) & is.na(age_at_death))),
    nrow(df %>% filter(language == lang, Gender == "Man",   !is.na(Gender) & is.na(age_at_death))),
    nrow(df %>% filter(language == lang, is.na(Gender) &  !is.na(age_at_death))),
    nrow(df %>% filter(language == lang, is.na(Gender) &   is.na(age_at_death)))
  )
}

# --- 3. Compute all groups ---
qc_fr    <- count_subsets(obituaries_quebec,  "fr")
qc_en    <- count_subsets(obituaries_quebec,  "eng")
ont_fr   <- count_subsets(obituaries_ontario, "fr")
ont_en   <- count_subsets(obituaries_ontario, "eng")

# --- 4. Combine results into table ---
NAtable <- cbind(
  Quebec_French  = unlist(qc_fr),
  Quebec_English = unlist(qc_en),
  Ontario_French = unlist(ont_fr),
  Ontario_English= unlist(ont_en)
)

# --- 5. Assign row/column names ---
rownames(NAtable) <- c(
  "Webscraped non-duplicates",
  "Age & Gender available",
  "Age & Gender - Women",
  "Age & Gender - Men",
  "Gender only (no age)",
  "Gender only - Women",
  "Gender only - Men",
  "Age only (no gender)",
  "Neither Age nor Gender"
)

colnames(NAtable) <- c("French", "English", "French", "English")

# --- Generate LaTeX table ---
xtable(NAtable,
       caption = "Extraction of demographic variables in French and English obituaries, by province, 2017–2022")

# The rest of the formatting was done in LaTeX



###############################################################################
# Table 4
###############################################################################

# --- Helper function: get total, retained, and % retention ---
get_counts <- function(df, lang) {
  total <- df %>%
    filter(language == lang) %>%
    group_by(Year) %>%
    summarise(total = n(), .groups = "drop")
  
  retained <- df %>%
    filter(language == lang, !is.na(Gender), !is.na(age_at_death)) %>%
    group_by(Year) %>%
    summarise(retained = n(), .groups = "drop")
  
  full <- left_join(total, retained, by = "Year") %>%
    mutate(
      retention_rate = round(100 * retained / total, 1)
    )
  
  return(full)
}

# --- Collect results for all 4 groups ---
qc_fr  <- get_counts(obituaries_quebec,  "fr")
qc_en  <- get_counts(obituaries_quebec,  "eng")
ont_fr <- get_counts(obituaries_ontario, "fr")
ont_en <- get_counts(obituaries_ontario, "eng")

# --- Combine into one table ---
tableyear <- data.frame(
  Year = 2017:2022,
  QC_French_Total     = qc_fr$total,
  QC_French_Retention = qc_fr$retention_rate,
  QC_English_Total     = qc_en$total,
  QC_English_Retention = qc_en$retention_rate,
  ONT_French_Total     = ont_fr$total,
  ONT_French_Retention = ont_fr$retention_rate,
  ONT_English_Total     = ont_en$total,
  ONT_English_Retention = ont_en$retention_rate
)

# --- Rename columns for clarity ---
colnames(tableyear) <- c("Year",
                         "QC-FR Total", "QC-FR % Retained",
                         "QC-EN Total", "QC-EN % Retained",
                         "ON-FR Total", "ON-FR % Retained",
                         "ON-EN Total", "ON-EN % Retained")

# --- Generate LaTeX table ---
xtable(tableyear,
       caption = "Retention rates (%) after demographic extraction by language and province, 2017–2022")

# The rest of the formatting was done in LaTeX



###############################################################################
# Table 5-6
###############################################################################

# --- Reset environment but keep chmd_deaths and chmd_exposures ---
rm(list = setdiff(ls(), c("chmd_deaths", "chmd_exposures", "chmd_life_tables",
                          "obits_life_tables")))

# --- 1. Load Quebec Obituaries Data ---
load("./Data/Quebec_obituaries_extracted.RData")


# --- 2. Mean age at death ---

# Mean age from obituaries
mean_age_obit <- obituaries_quebec %>% 
  group_by(Year, Gender) %>% 
  summarise(mean_age_obituaries = mean(age_at_death, na.rm = TRUE))

# Mean age from CHMD (weighted)
mean_age_CHMD <- chmd_deaths %>% 
  group_by(Year, Gender) %>% 
  summarise(mean_age_CHMD = weighted.mean(Age, Deaths))

# --- 3. Median age at death ---

# From obituaries
median_obituaries <- obituaries_quebec %>% 
  group_by(Year, Gender) %>% 
  summarise(median_obituaries = median(age_at_death, na.rm = TRUE))


median_CHMD <- chmd_deaths %>%
  select(-Total) %>%
  mutate(Deaths = round(Deaths)) %>%
  uncount(Deaths) %>%
  group_by(Year, Gender) %>%
  summarise(median_CHMD = median(Age))

# --- 4. Merge and calculate differences ---

median_age <- left_join(median_obituaries, median_CHMD, by = c("Year", "Gender")) %>%
  mutate(difference = median_obituaries - median_CHMD)

mean_age <- left_join(mean_age_obit, mean_age_CHMD, by = c("Year", "Gender")) %>%
  mutate(difference = mean_age_obituaries - mean_age_CHMD)

# --- 5. Create separate tables for Women and Men ---

# Table for Women
wtable <- cbind(
  mean_age %>% filter(Gender == "Woman") %>% select(Year, mean_age_obituaries, mean_age_CHMD, difference),
  median_age %>% filter(Gender == "Woman") %>% ungroup() %>% select(median_obituaries, median_CHMD, difference)
)
colnames(wtable) <- c("Year", "Mean obituaries", "Mean CHMD", "delta1",
                      "Median obituaries", "Median CHMD", "delta2")

wtable_x <- xtable(wtable,
                   caption = "Mean and Median Age at Death Comparison: Women, Quebec, 2017–2022")

print(wtable_x,
      include.rownames = FALSE,
      caption.placement = "top")

# Table for Men
mtable <- cbind(
  mean_age %>% filter(Gender == "Man") %>% select(Year, mean_age_obituaries, mean_age_CHMD, difference),
  median_age %>% filter(Gender == "Man") %>% ungroup() %>% select(median_obituaries, median_CHMD, difference)
)
colnames(mtable) <- c("Year", "Mean obituaries", "Mean CHMD", "delta1",
                      "Median obituaries", "Median CHMD", "delta2")

mtable_x <- xtable(mtable,
                   caption = "Mean and Median Age at Death Comparison: Men, Quebec, 2017–2022")

print(mtable_x,
      include.rownames = FALSE,
      caption.placement = "top")


###############################################################################
# Tables 7–8: Modal Age at Death and Life Expectancy Comparisons
###############################################################################

# --- 1. Merge CHMD and Obituary Data: Adjust obituary deaths to CHMD totals ---

df_combined <- left_join(
  chmd_deaths %>% 
    group_by(Gender, Year) %>% 
    mutate(Total = sum(Deaths), pct = Deaths / Total) %>%
    select(-Deaths),
  obituaries_quebec %>% 
    group_by(Gender, Year, age_at_death) %>% 
    summarise(Deaths = n(), .groups = "drop") %>% 
    mutate(Total = sum(Deaths, na.rm = TRUE))
) %>%
  mutate(Deaths = Total * pct)

df_combined <- left_join(
  chmd_deaths %>% 
    group_by(Gender, Year) %>% 
    mutate(Total = sum(Deaths)),
  obituaries_quebec %>% 
    group_by(Gender, Year, age_at_death) %>% 
    summarise(Deaths_obit = n(), .groups = "drop") %>%
    mutate(Total = sum(Deaths_obit), pct = Deaths_obit / Total) %>%
    rename(Age = age_at_death) %>%
    select(Age, Gender, Year, pct)
) %>%
  mutate(
    pct = ifelse(is.na(pct), 0, pct),
    Deaths_obit = Total * pct
  )

# --- 2. Attach Exposure Data ---
df_combined <- left_join(df_combined, chmd_exposures)

# --- 3. Reshape and label data sources ---
df_combined <- df_combined %>%
  pivot_longer(c("Deaths", "Deaths_obit"), names_to = "Source", values_to = "Deaths") %>%
  mutate(Source = ifelse(Source == "Deaths", "CHMD", "Obituaries"))

# --- 4. Calculate modal age at death using smoothed distribution ---
library(MortalitySmooth)

M_function <- function(source, year, gender) {
  
  subset_df <- df_combined %>% 
    filter(Gender == gender, Year == year, Age >= 5, Source == source)
  
  D <- subset_df$Deaths
  E <- subset_df$Exposures
  ages <- subset_df$Age
  
  fit <- Mort1Dsmooth(x = ages + 0.5, y = D, offset = log(E))
  
  delta <- 0.01
  age.narrow <- seq(from = min(ages), to = max(ages) + 1 - delta, by = delta)
  
  log.mu.hat <- predict(fit, newdata = age.narrow)
  mu.hat <- exp(log.mu.hat)
  l.hat <- exp(-cumsum(mu.hat * delta))
  d.hat <- mu.hat * l.hat
  
  return(age.narrow[which.max(d.hat)])
}

# --- 5. Calculate modal age for Men and Women ---

years <- 2017:2022

M_males <- data.frame(
  Year = years,
  M_obit = sapply(years, function(y) M_function("Obituaries", y, "Man")),
  M_stats = sapply(years, function(y) M_function("CHMD", y, "Man"))
)

M_females <- data.frame(
  Year = years,
  M_obit = sapply(years, function(y) M_function("Obituaries", y, "Woman")),
  M_stats = sapply(years, function(y) M_function("CHMD", y, "Woman"))
)

# Combine modal age data
M_s <- bind_rows(
  M_females %>% mutate(Gender = "Woman"),
  M_males %>% mutate(Gender = "Man")
) %>%
  mutate(diff_mode = M_obit - M_stats)

# --- 6. Merge with life expectancy data ---

le_diff_df <- left_join(
  chmd_life_tables %>% select(Age, Year, Gender, e_x),
  obits_life_tables %>% select(Age, Year, Gender, e_x) %>% rename(e_x_obit = e_x),
  by = c("Age", "Year", "Gender")
) %>%
  mutate(diff = e_x_obit - e_x)

# Final dataset at age 0 (life expectancy at birth)
final_combined <- left_join(
  M_s,
  le_diff_df %>% filter(Age == 0) %>% select(-Age),
  by = c("Year", "Gender")
) %>%
  select(Year, Gender, M_stats, M_obit, diff_mode, e_x, e_x_obit, diff)

# --- 7. Format final tables with xtable ---
library(xtable)

# Combined table (Men and Women)
final_table <- xtable(
  final_combined,
  caption = "Life Expectancy at Birth and Modal Age at Death Differences by Data Source and Gender, Quebec, 2017–2022"
)

# Women
wtable <- final_combined %>%
  filter(Gender == "Woman") %>%
  select(Year, e_x, e_x_obit, diff, M_stats, M_obit, diff_mode)

wtable_x <- xtable(
  wtable,
  caption = "Life Expectancy at Birth and Modal Age at Death Differences by Data Source, Women, Quebec, 2017–2022"
)

print(wtable_x,
      caption.placement = "top",
      include.rownames = FALSE)

# Men
mtable <- final_combined %>%
  filter(Gender == "Man") %>%
  select(Year, e_x, e_x_obit, diff, M_stats, M_obit, diff_mode)

mtable_x <- xtable(
  mtable,
  caption = "Life Expectancy at Birth and Modal Age at Death Differences by Data Source, Men, Quebec, 2017–2022"
)

print(mtable_x,
      caption.placement = "top",
      include.rownames = FALSE)



# For any questions please e-mail
# pietro.violo@umontreal.ca