library(tidyverse)
library(WikidataQueryServiceR)
library(lubridate)
library(stringr)
library(magrittr)
library(plotly)
library(DT)
library(scales)
library(ggthemes)


my_wiki_query <- query_wikidata('SELECT ?item ?dob ?dod ?sex_or_gender WHERE {
  ?item wdt:P31 wd:Q5;
  OPTIONAL { ?item wdt:P570 ?dod }
  OPTIONAL { ?item wdt:P569 ?dob }
  OPTIONAL { ?item wdt:P21 ?sex_or_gender. }
}')



my_wiki_query <- readRDS("./data/my_wiki_query.rds")

count_query <- nrow(my_wiki_query)


# sex notation
my_wiki_query %<>% 
  mutate(simplified_sex = case_when(
    sex_or_gender == "http://www.wikidata.org/entity/Q6581072" ~ 1, #female
    sex_or_gender == "http://www.wikidata.org/entity/Q6581097" ~ 2, #male
    TRUE ~ 3 #not female, not male
    ))

my_wiki_query %<>% 
  dplyr::distinct(item, .keep_all= TRUE) %>% 
  dplyr::mutate(full_year_of_birth = str_extract(dob, "-?[:digit:]{4}")) %>% 
  dplyr::mutate(full_year_of_death = str_extract(dod, "-?[:digit:]{4}")) %>% 
  dplyr::mutate(approx_age_at_death = as.numeric(full_year_of_death) - as.numeric(full_year_of_birth)) 

my_wiki_query$full_year_of_birth <- as.numeric(my_wiki_query$full_year_of_birth) 
my_wiki_query$full_year_of_death <- as.numeric(my_wiki_query$full_year_of_death)   
my_wiki_query$simplified_sex <- as_factor(my_wiki_query$simplified_sex)

my_wiki_query %<>% mutate(decade_birth = full_year_of_birth %/% 10)


my_wiki_query %>% 
  filter(full_year_of_birth >= 1000 & full_year_of_birth < 2022) %>% 
  filter(simplified_sex == 1 | simplified_sex == 2) %>% 
  ggplot(aes(x = full_year_of_birth)) +
  geom_histogram(binwidth = 10) +
  scale_y_log10(breaks=c(100, 1000, 10000, 100000), labels = c("100", "1000", "10,000", "100,000")) +
  scale_x_continuous(breaks = seq(1000, 2021, 100)) +
  theme_economist() + 
  labs(title = "Distribution of births per year",
       x = "Year",
       y = "Birth count per year",
       subtitle = "Years 1000 to 2021\nLog Scale Y-axis", 
       caption = "Including only records with birth years,\nand restricted to sex of either male or female.") +
  coord_cartesian(ylim = c(100, 200000))


male_color <- "#2E45B8"
female_color <- "#C91D42"

sex_proportion_plot <- my_wiki_query %>% 
  filter(simplified_sex == 1 | simplified_sex == 2) %>% 
  filter(full_year_of_birth >= 1000 & full_year_of_birth < 2021) %>% 
#  group_by(decade_birth) %>% 
  ggplot(aes(x = full_year_of_birth, fill = simplified_sex, color = simplified_sex)) +
  geom_bar(position = "fill") +
  scale_fill_manual(labels = c("Female", "Male"), values = c(female_color, male_color)) +
  scale_color_manual(labels = c("Female", "Male"), values = c(female_color, male_color)) +
  labs(title = "Proportion of men and women over time in birth cohorts", 
       x = "Year of birth",
       y = " ", 
       fill = "Sex", 
       color = NULL) +
  geom_hline(yintercept = 0.5) +
  theme_economist() +
  scale_y_continuous(breaks=c(0, 0.25, .5, 0.75, 1), 
                     labels = c(" ", " ", "50%", " ", " ")) +
  scale_x_continuous(breaks = seq(1000, 2021, 100)) +
  guides(color = "none")

sex_proportion_plot


ce_dates <- my_wiki_query %>% 
  filter(full_year_of_birth >= 1000 & full_year_of_birth <= 2021) %>% 
  filter(full_year_of_death >= 1000 & full_year_of_birth <= 2021) %>% 
  mutate(date_of_birth = as_date(dob)) %>%
  mutate(year_birth = year(date_of_birth)) %>%
  mutate(month_birth = month(date_of_birth)) %>%
  mutate(day_birth = day(date_of_birth)) %>%
  mutate(date_of_death = as_date(dod)) %>%
  mutate(year_death = year(date_of_death)) %>%
  mutate(month_death = month(date_of_death)) %>%
  mutate(day_death = day(date_of_death)) %>%
  mutate(lifespan = as.duration(interval(date_of_birth, date_of_death))) %>% 
  filter(lifespan < years(123)) %>% 
  filter(lifespan >= years(0))


average_lifespan_by_sex <- ce_dates %>% 
  group_by(year_birth, simplified_sex) %>% 
  summarise(average_age = mean((lifespan)))

average_lifespan <- ce_dates %>% 
  group_by(year_birth) %>% 
  summarise(average_age = mean((lifespan)))

plot_average_lifespan <- average_lifespan %>%  
  ggplot(aes(x = year_birth, y = (average_age/31557600))) +
  geom_smooth(color = 	"#F97A1F") +
  geom_point(alpha = 0.05, color = 	"#F97A1F") +
  labs(title = "Average age of death over time",
       x = "Birth Year Cohort",
       y = "Age") 

plot_average_lifespan + theme_economist() + scale_x_continuous(breaks = seq(1000, 2021, 100))

average_lifespan_by_sex <- ce_dates %>% 
  group_by(year_birth, simplified_sex) %>% 
  summarise(average_age = mean((lifespan)))

plot_average_lifespan_by_sex <- average_lifespan_by_sex %>% 
  filter(simplified_sex != 3) %>% 
  ggplot(aes(x = year_birth, y = (average_age/31557600), color = simplified_sex)) +
  geom_smooth() +
  geom_point(alpha = 0.05) +
  scale_color_manual(labels = c("Female", "Male"), values = c(female_color, male_color)) +
  labs(title = "Average age of death over time by sex",
       x = "Birth Year Cohort",
       y = "Age", 
       color = "Sex") 

plot_average_lifespan_by_sex + theme_economist() + scale_x_continuous(breaks = seq(1000, 2021, 100)) 


average_lifespan_by_sex_dc <- ce_dates %>% 
  group_by(year_death, simplified_sex) %>% 
  summarise(average_age = mean((lifespan)))

average_lifespan_dc <- ce_dates %>% 
  group_by(year_death) %>% 
  summarise(average_age = mean((lifespan)))

plot_average_lifespan_dc <- average_lifespan_dc %>%  
  ggplot(aes(x = year_death, y = (average_age/31557600))) +
  geom_smooth(color = 	"#F97A1F") +
  geom_point(alpha = 0.05, color = 	"#F97A1F") +
  labs(title = "Average age of death over time",
       x = "Death Year Cohort",
       y = "Age") 

plot_average_lifespan_dc  + theme_economist() + scale_x_continuous(breaks = seq(1000, 2021, 100))

plot_average_lifespan_by_sex_dc <- average_lifespan_by_sex_dc %>% 
  filter(simplified_sex != 3) %>% 
  ggplot(aes(x = year_death, y = (average_age/31557600), color = simplified_sex)) +
  geom_smooth() +
  geom_point(alpha = 0.05) +
  scale_color_manual(labels = c("Female", "Male"), values = c(female_color, male_color)) +
  labs(title = "Average age of death over time by sex",
       x = "Death Year Cohort",
       y = "Age", 
       color = "Sex") 

plot_average_lifespan_by_sex_dc + theme_economist() + scale_x_continuous(breaks = seq(1000, 2021, 100))


average_lifespan_dc_20 <- ce_dates %>% 
  filter(lifespan > 630720000) %>% 
  group_by(year_death) %>% 
  summarise(average_age = mean((lifespan)))

plot_average_lifespan_dc_20 <- average_lifespan_dc_20  %>%  
  ggplot(aes(x = year_death, y = (average_age/31557600))) +
  geom_smooth(color = 	"#F97A1F") +
  geom_point(alpha = 0.05, color = 	"#F97A1F") +
  labs(title = "Average age of death over time",
       subtitle = "Only lifespans greater than 20 years",
       x = "Death Year Cohort",
       y = "Age") 

plot_average_lifespan_dc_20 + theme_economist() + scale_x_continuous(breaks = seq(1000, 2021, 100))


plot_average_lifespan_dc <- average_lifespan_dc %>%  
  filter(year_death >= 1600) %>% 
  ggplot(aes(x = year_death, y = (average_age/31557600))) +
  geom_smooth(color = 	"#F97A1F") +
  geom_point(alpha = 0.05, color = 	"#F97A1F") +
  labs(title = "Average age of death over time",
       x = "Death Year Cohort",
       y = "Age") 

plot_average_lifespan_dc + 
  theme_economist() + 
  scale_x_continuous(breaks = seq(1000, 2021, 100)) + 
    coord_cartesian(ylim = c(40, 90))

plot_average_lifespan_by_sex_dc <- average_lifespan_by_sex_dc %>% 
  filter(year_death >= 1600) %>% 
  filter(simplified_sex != 3) %>% 
  ggplot(aes(x = year_death, y = (average_age/31557600), color = simplified_sex)) +
  geom_smooth() +
  geom_point(alpha = 0.05) +
  scale_color_manual(labels = c("Female", "Male"), values = c(female_color, male_color)) +
  labs(title = "Average age of death over time",
       x = "Death Year Cohort",
       y = "Age", 
       color = "Sex") +
  guides(color = "none")

plot_average_lifespan_by_sex_dc + 
  theme_economist() + 
  scale_x_continuous(breaks = seq(1000, 2021, 100)) + 
  coord_cartesian(ylim = c(40, 90))


year_cutoff <- 1800

plot_average_lifespan_dc <- average_lifespan_dc %>%  
  filter(year_death >= year_cutoff) %>% 
  ggplot(aes(x = year_death, y = (average_age/31557600))) +
  geom_line(color = 	"#F97A1F") +
 # geom_point(alpha = 0.05, color = 	"#F97A1F") +
  labs(title = "Average age of death over time",
       x = "Death Year Cohort",
       y = "Age") 

plot_average_lifespan_dc + 
  theme_economist() + 
  scale_x_continuous(breaks = seq(1000, 2021, 100)) + 
    coord_cartesian(ylim = c(40, 90))

plot_average_lifespan_by_sex_dc <- average_lifespan_by_sex_dc %>% 
  filter(year_death >= year_cutoff) %>% 
  filter(simplified_sex != 3) %>% 
  ggplot(aes(x = year_death, y = (average_age/31557600), color = simplified_sex)) +
  geom_line() +
  scale_color_manual(labels = c("Female", "Male"), values = c(female_color, male_color)) +
  labs(title = "Average age of death over time by sex",
       x = "Death Year Cohort",
       y = "Age", 
       color = "Sex") +
  guides(color = "none")

plot_average_lifespan_by_sex_dc + 
  theme_economist() + 
  scale_x_continuous(breaks = seq(1000, 2021, 100)) + 
  coord_cartesian(ylim = c(40, 90))




heatmap_dates <- ce_dates %>% 
  mutate(pad_month_birth = str_pad(month_birth, width = 2, side = "left", pad = 0)) %>% 
  mutate(pad_day_birth = str_pad(day_birth, width = 2, side = "left", pad = 0)) %>% 
  unite("heatdate", c(year_birth, pad_month_birth, pad_day_birth), sep = "-", remove = FALSE) %>% 
  mutate(doy = lubridate::yday(as_date(heatdate)))

heatmap_dates_subset <- heatmap_dates %>% dplyr::select(year_birth, heatdate, doy)

heatmap_data <- heatmap_dates_subset %>% 
  group_by(year_birth, doy) %>% 
  summarise(count = n()) %>% 
  ungroup() %>% 
  group_by(year_birth) %>% 
  mutate(year_total = sum(count)) %>% 
  mutate(proportion_year = count/year_total) %>% 
  mutate(too_high = case_when(
    proportion_year > 0.003 ~ 1, #roughly 10% higher than expected
    TRUE ~ 0
  ))


too_many_births_plot <- heatmap_data  %>% 
  filter(year_birth >= 1500) %>% 
  ggplot(aes(x = year_birth, y = doy, fill = too_high, text = )) + 
  geom_tile() +
  scale_fill_gradient(low = "#475ED1", high = "#F9C31F") + # mid 
  theme_economist() +
  theme(legend.position = "none") +
  labs(x = "Year", 
       y = " ", 
       title = "Days of the year with more than expected birthdates",
       subtitle = "Yellow represents with \"too\" many births \nY-axis is the day of the year with tick marks to show the start of the month") +
  scale_y_continuous(breaks=c(1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335),
                     labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))

too_many_births_plot


too_many_births_plot <- heatmap_data  %>% 
  filter(doy < 31) %>% 
  ggplot(aes(x = year_birth, y = doy, fill = too_high, text = )) + 
  geom_tile() +
  scale_fill_gradient(low = "#475ED1", high = "#F9C31F") + # mid 
  theme_economist() +
  theme(legend.position = "none") +
  labs(x = "Year", 
       y = " ", 
       title = "Days of the year with more than expected birthdates",
       subtitle = "Yellow represents with \"too\" many births \nY-axis is the day of the year with tick marks to show the start of the month") +
  scale_y_continuous(breaks=c(1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335),
                     labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))

too_many_births_plot




sessionInfo()


image_for_blog_post <- plot_average_lifespan_by_sex_dc + 
    theme_economist() + 
     scale_x_continuous(breaks = seq(1000, 2021, 100)) + 
     coord_cartesian(ylim = c(40, 90))

ggsave(file = "./images/image_for_blog_post.png" )