lpath =  paste("C://Users/", Sys.info(  )["login"], "/R.lib", sep = "")
.libPaths(lpath)

library(tidyverse)
## -- Attaching packages ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.2.1     v purrr   0.3.3
## v tibble  2.1.3     v dplyr   0.8.4
## v tidyr   1.0.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0
## -- Conflicts ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(knitr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(here)
## here() starts at C:/Users/LayC/Projects/cira3-mortality
## 
## Attaching package: 'here'
## The following object is masked from 'package:lubridate':
## 
##     here
library(readxl)
library(mvmeta)
## This is mvmeta 1.0.3. For an overview type: help('mvmeta-package').
library(dlnm)
## This is dlnm 2.3.9. For details: help(dlnm) and vignette('dlnmOverview').
if(!dir.exists(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"))){ dir.create(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), recursive = TRUE)}

Goals

  1. To check the MCMC output for projected - hindcast AN from project_spline_mcmc.r called by CIRA-3_Mortality_-Base_y4_load-precalc_CLay.Rmd The function takes a very time to run, as calculating projected-hindcast for each simulated set of coefficients and each individual day in the projected period provides the most flexibility to calculate the different windows (seasonal, temperature range, day of year).
  2. Compare the detailed single city output intended for delivery and example data, values from the MCMC process and, the overall central estimate (proj-hind) values used in mapping to ensure they match.
  3. Create plots and tables for the manuscript
load(here("output/y4-final-manuscript-output/an_ci_dtC2_yr_season.Rdata")) # MCMC based estimates, with CI, for seasons, years, and temperature ranges; dtC2 and dtC 1 are still running

full_pop2010 <- read_xls(here("data/City-County-Cluster-2014-06-24_POP2010.xls"), "City-County-Cluster-2014-06-24_") # 2010 populations used for normalization
raw_boston_base <- read_xlsx(here("output/y4-final-manuscript-output/2019-02-05_projected_mortality_2003-2013-dr.xlsx"), "baseline") # baseline central estimate  
raw_boston_dtc6 <- read_xlsx(here("output/y4-final-manuscript-output/2019-02-05_projected_mortality_2003-2013-dr.xlsx"), "dtC_6") # dt_C central estimate  
summary_out_mean <- read_csv(here("output/y4-final-manuscript-output/2020-02-04_city_dT_AN_y4_dif_proj-hind_summary.csv")) # summary central estimate file
## Parsed with column specification:
## cols(
##   cluster = col_double(),
##   city = col_character(),
##   CITYNAME = col_character(),
##   FIPS_txt = col_character(),
##   dT_c = col_double(),
##   baseline_proj_mort = col_double(),
##   sum_an_total = col_double(),
##   sum_an_total_bnd = col_double()
## )

Ensure all cities present full_pop2010

full_pop2010 %>% mutate(CITYNAME = casefold(CITYNAME) , FIPStext = FIPStext) %>% 
  select(CITYNAME, City) %>%  unique () %>% 
  anti_join(city_ci_sums [, c("CITYNAME", "FIPS_txt" )])
## Joining, by = "CITYNAME"
full_pop2010 <- full_pop2010 %>% mutate(CITYNAME = casefold(CITYNAME) , FIPS_txt = FIPStext) 
city_ci_sums [, c("CITYNAME", "FIPS_txt" )] %>% anti_join(full_pop2010) %>% unique()
## Joining, by = c("CITYNAME", "FIPS_txt")
## all good
city_ci_sums %>% filter(!is.na(CITYNAME)) %>% pull(CITYNAME) %>% unique() %>% length() # 208 = yes
## [1] 208

Goal 1

Check the MCMC output for match against the values used in mapping. First normalize by 2010 population.

city_ci_sums_pop <-  full_pop2010 %>% 
  mutate(pop_over65 = USA_Counties.AGE_65_74 + USA_Counties.AGE_75_84 + USA_Counties.AGE_85_UP) %>% 
  group_by(CITYNAME, Cluster ) %>% 
  summarise(n_counties = n() , 
            pop2010 = sum(USA_Counties.POP2010), 
            pop_over65 = sum(pop_over65)) %>%
  right_join(city_ci_sums) %>%
  filter(!is.na(CITYNAME)) %>%  
  pivot_longer(cols = mean_val:last_col()) %>% 
  mutate(value = as.numeric(value)) %>%
  mutate (an_pop2010_100k = case_when(grepl("an_", tm_var)  ~ (value/pop2010)*100000,
                                   TRUE ~ value) ) %>%
  ungroup()
## Joining, by = "CITYNAME"
nrow(city_ci_sums) 
## [1] 222193
city_ci_sums_pop  

Goal 2: Check Boston delivery file

Check calculations for consistency with both MCMC and proj-hind used in mapping

boston_base_yr_mean <- raw_boston_base %>% 
  select(city, CITYNAME, cluster, year, gcm, dT_c  , FIPS_txt, doy, tmean0, tmean15, rr_tmean0:rr_total,
         tmean0_an_2010_100k:tot_an_2010_100k, an_tmean0:an_total ) %>%
  pivot_longer(-(city:doy)) %>% 
  mutate( type = case_when (grepl( "an_", name ) ~ "mean_year_sum"  , TRUE ~ "mean_year_mean" )) %>% 
  group_by(city, CITYNAME, cluster, year, gcm, dT_c  , FIPS_txt, type, name) %>% 
  summarise(yr_sum =  sum(value) , yr_mean =mean(value) ) %>% 
  ungroup()%>% 
  mutate(value = case_when (grepl( "an_", name ) ~ yr_sum
                                       , TRUE ~ yr_mean))%>%
  group_by(city, CITYNAME, cluster,  FIPS_txt, name, type) %>%
  summarise(value = mean(value))%>% ungroup()


boston_dtc6_yr_mean <- raw_boston_dtc6 %>% 
  select(city, CITYNAME, cluster, year, gcm, dT_c  , FIPS_txt, doy, tmean0, tmean15, rr_tmean0:rr_total,tmean0_an_2010_100k:tot_an_2010_100k, an_tmean0:an_total )%>%
  pivot_longer(-(city:doy)) %>% 
  mutate( type = case_when (grepl( "an_", name ) ~ "mean_year_sum"  , TRUE ~ "mean_year_mean" )) %>% 
  group_by(city, CITYNAME, cluster, year, gcm, dT_c  , FIPS_txt, type, name) %>% 
  summarise(yr_sum =  sum(value) , yr_mean =mean(value) ) %>% 
  ungroup()%>% 
  mutate(value = case_when (grepl( "an_", name ) ~ yr_sum
                                       , TRUE ~ yr_mean))%>%
  group_by(city, CITYNAME, cluster,  FIPS_txt, name, type) %>%
  summarise(value = mean(value))%>% ungroup()

boston_base_yr_mean  %>% 
  left_join(boston_dtc6_yr_mean, by = c("city", "CITYNAME", "cluster", "FIPS_txt", "name", "type" ), suffix = c("_base", "_dt6") ) %>% 
  mutate(value_dif = value_dt6 - value_base)%>% ### the raw sumamry matches that in MATCHES 2020-02-04-city_dT_AN_y4_dif_proj-hind_summary.csv (whole number match; rounding off)
  filter(name %in% c("an_total" )) 
summary_out_mean %>% filter(CITYNAME =="boston")

Goal 3

Create data for plotting and tables. Divide the MCMC output into the aggregates for year, season, and temperature range.

city_ci_sums_year <- city_ci_sums_pop %>% 
  filter(!is.na(dT_c) & season ==12) %>%  
  mutate( Cluster= factor(cluster, levels = 1:9,  labels = c("1: NE Coast", "2: NE", "3: Midwest", "4: SE Central", "5: W Coast", 
                                                                       "6: N Gulf", "7: SE Reach", "8: SW", "9: West" )), 
          `Historical period` = factor(time_period, levels = c(1,2,3,4), 
                                       labels = c("1973-1983", "1983-1993", "1993-2003", "2003-2013") ) )

city_ci_sums_season <- city_ci_sums_pop %>% filter(!is.na(dT_c) & season %in% 1:4) %>% 
  mutate( Cluster= factor(cluster, levels = 1:9,  labels = c("1: NE Coast", "2: NE", "3: Midwest", "4: SE Central", "5: W Coast", 
                                                                       "6: N Gulf", "7: SE Reach", "8: SW", "9: West" )), 
          `Historical period` = factor(time_period, levels = c(1,2,3,4), 
                                       labels = c("1973 -1983", "1983 -1993", "1993 -2003", "2003 -2013") ),
          Season = factor(season, levels = c(1:4),
                         labels = c("Winter\n(Dec-Feb)", "Spring\n(Mar-June)", "Summer\n(July-Sep)", "Fall\n(Oct-Nov)")))

city_ci_sums_tmprng <- city_ci_sums_pop %>% filter(!is.na(dT_c) & season == 0) %>% 
   mutate( Cluster= factor(cluster, levels = 1:9,  labels = c("1: NE Coast", "2: NE", "3: Midwest", "4: SE Central", "5: W Coast", 
                                                                       "6: N Gulf", "7: SE Reach", "8: SW", "9: West" )), 
          `Historical period` = factor(time_period, levels = c(1,2,3,4), 
                                       labels = c("1973-1983", "1983-1993", "1993-2003", "2003-2013") ),
          `Range C` = factor(tmp_rng, levels = c(-1, -0.5,0 , 0.5, 1), 
                             labels = c("Coldest", "Cold", "Mid", 
                                        "Hot","Hottest" ))) 

Goal 2: Check consistency between MCMC central estimate and central estimate in proj-hind

This should be very similar, within rounding error. If it is not, it indicates an error somewhere in the aggregation of the mean value for AN.

difs <- city_ci_sums_year %>% select(-Cluster,  -FIPS_txt) %>% 
  filter(  name == "mean_val" & tm_var == "an_total" & time_period == 4) %>% 
  left_join(summary_out_mean) %>% 
  select(cluster, CITYNAME ,city,  dT_c, pop2010, tm_var, value,sum_an_total, n_samp) %>% 
  mutate( mcmc_mean_an_2010_norm =( value/pop2010)*100000, 
          central_est_an_2010_norm = (sum_an_total/pop2010)*100000 , # normalize by 2010 population
         dif_2010_norm =  central_est_an_2010_norm-mcmc_mean_an_2010_norm  ) # get difference
## Joining, by = c("CITYNAME", "dT_c", "cluster")
difs
difs%>% # slight differences.  Nothing that greatly affects output.
  ggplot(aes(dif_2010_norm)) + 
  geom_histogram() +
  xlab("Difference\nMCMC mean AN - Central estimate AN (per 100k 2010 pop)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

difs%>%
  ggplot(aes(x = cluster, y = dif_2010_norm, color = dT_c )) + geom_point()  +
  ylab("Difference\nMCMC mean AN - Central estimate AN (per 100k 2010 pop)")

Goal 3

Create table for manuscript, per Marcus Sarofim request for a table of national totals with lower and upper range. Note that this is normalized to 2010 population.

tp_adapt_sum <- city_ci_sums_year %>% filter(tm_var == "an_total" ) %>% 
  select(Cluster, CITYNAME, FIPS_txt, `Historical period`, `Degrees C change` = dT_c, tm_var, n_samp, name, an_pop2010_100k) %>%
  pivot_wider(names_from = "name", values_from = "an_pop2010_100k") %>%
  group_by(`Degrees C change`, `Historical period`  ) %>%
    summarise(`National change in deaths due to temperature ` = round(sum(mean_val),0),
              `Lower Range` =  round(sum(lci_66),0),
              `Upper Range` =  round(sum(uci_66),0))  
tp_adapt_sum %>% 
  write.csv(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/summary_AN_tp-adapt_dif.csv"))
tp_adapt_sum %>% knitr::kable( padding = 2)
Degrees C change Historical period National change in deaths due to temperature Lower Range Upper Range
2 1973-1983 1352 553 2097
2 1983-1993 456 -42 919
2 1993-2003 29 -355 389
2 2003-2013 -129 -474 217
3 1973-1983 2669 1433 3833
3 1983-1993 990 197 1725
3 1993-2003 267 -320 822
3 2003-2013 -71 -645 492
4 1973-1983 4282 2492 6056
4 1983-1993 1719 642 2787
4 1993-2003 690 -147 1499
4 2003-2013 114 -723 941
5 1973-1983 6389 3766 8920
5 1983-1993 2744 1152 4337
5 1993-2003 1359 78 2594
5 2003-2013 451 -755 1660
6 1973-1983 10402 7860 12989
6 1983-1993 4825 2953 6665
6 1993-2003 2883 1428 4331
6 2003-2013 1280 -371 2944
tp_y4_sum <-  city_ci_sums_year %>% filter(tm_var == "an_total"& time_period == 4 ) %>% 
  select(Cluster, CITYNAME, FIPS_txt, `Degrees C change` = dT_c, tm_var, n_samp, name, an_pop2010_100k) %>%
  pivot_wider(names_from = "name", values_from = "an_pop2010_100k") %>%
  group_by(`Degrees C change`  ) %>%
    summarise(`National change in deaths due to temperature ` = round(sum(mean_val), 0),
              `Lower Range` = round(sum(lci_66), 0),
              `Upper Range` = round(sum(uci_66), 0)) 
tp_y4_sum%>% 
  write.csv(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/summary_AN_y4_dif.csv"))
tp_y4_sum %>% knitr::kable( padding = 2)
Degrees C change National change in deaths due to temperature Lower Range Upper Range
2 -129 -474 217
3 -71 -645 492
4 114 -723 941
5 451 -755 1660
6 1280 -371 2944

Make Plots of MCMC results

First, plot attributatble risk for projection - hindcast relative to season and the time period of the model fit. This shows adaptation over time, and how it differs by cluster.

for (i in 2:6) {# i == 2
ci_season_rr_2010100k_tp_facet <- city_ci_sums_season %>% 
  filter(season %in% c(1:4) & grepl("rr_total", tm_var) & dT_c == i) %>% # select out just the seasonal values and rr
  select( -an_pop2010_100k  ) %>% 
  pivot_wider() %>% 
  ggplot(aes(x = Season, 
             y = mean_val, 
             color = Cluster, 
             group = Cluster)) +
  geom_point(   position = position_dodge(width = 0.8), size = 0.9 ) + 
  geom_errorbar(aes(x=Season, ymin = lci_66, ymax = uci_66), width = 0  ,  
                position = position_dodge(width = 0.8) )+ 
  scale_color_viridis_d(option = "plasma" ) + 
  ylab(bquote('Seasonal AR'[Delta*degree*C==.(i)]))+  # paste0("Total Attributable Risk\nRelative to Hindcast at ",i ," Δ°C" ))  
  theme_classic() +
  facet_wrap( vars(`Historical period`))

print(ci_season_rr_2010100k_tp_facet)
 
 
ggsave(  paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, "_rr_season_ci_2010100k_tp_facet.pdf" ), plot = ci_season_rr_2010100k_tp_facet )
ggsave( paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, "_rr_season_ci_2010100k_tp_facet.png"), plot = ci_season_rr_2010100k_tp_facet )}
## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

Next plot attributable risk relative to hindcast against year; this is to have an option for comparison.

for (i in 2:6) {# i = 2
  rr_season_ci_2010100k_seas_facet <- city_ci_sums_season %>% 
  filter(season %in% c(1:4) & grepl("rr_total", tm_var)& dT_c == i) %>% # select out just the seasonal values and rr
  select( -an_pop2010_100k  ) %>% 
  pivot_wider() %>% 
  #mutate(`Historical period` = str_wrap(`Historical period`, width = 4)) %>% 
  ggplot(aes(x = `Historical period`, 
             y = mean_val, 
             color = Cluster, 
             group = Cluster)) +
  geom_point(   position = position_dodge(width = 0.8) , size = 0.9) + 
  geom_errorbar(aes(x=`Historical period`, 
                    ymin = lci_66, ymax = uci_66), width = 0  ,  position = position_dodge(width = 0.8) )+ 
   ylab(bquote('Seasonal AR '[Delta*degree*C==.(i)])) +#   paste0("Total Attributable Risk\nRelative to Hindcast at  ",i ," Δ°C")) +
  scale_x_discrete(labels = function(`Historical period`) str_wrap(`Historical period` , width = 4) ) +
  scale_color_viridis_d(option = "plasma" ) + 
  theme_classic() +
  facet_wrap( vars(Season) )
 
print(rr_season_ci_2010100k_seas_facet)
 
 
ggsave( paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, "_rr_season_ci_2010100k_seas_facet.pdf" ), plot = rr_season_ci_2010100k_seas_facet )

ggsave( paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, "_rr_season_ci_2010100k_seas_facet.png"), plot = rr_season_ci_2010100k_seas_facet )}
## Saving 6 x 4 in image
## Saving 6 x 4 in image

## Saving 6 x 4 in image
## Saving 6 x 4 in image

## Saving 6 x 4 in image
## Saving 6 x 4 in image

## Saving 6 x 4 in image
## Saving 6 x 4 in image

## Saving 6 x 4 in image
## Saving 6 x 4 in image

Also plot the attributable mortality, which accounts for seasonal differences in mortality not related to temperature.

for (i in 2:6) {
  ci_season_rr_2010100k_tp_facet <- city_ci_sums_season %>% 
  filter(season %in% c(1:4) & grepl("an_total", tm_var) & dT_c == i) %>% # select out just the seasonal values and rr
  mutate(value = an_pop2010_100k) %>% select( -an_pop2010_100k  ) %>% 
  pivot_wider() %>% 
  ggplot(aes(x = Season, 
             y = mean_val, 
             color = Cluster, 
             group = Cluster)) +
  geom_point(   position = position_dodge(width = 0.8), size = 0.9 ) + 
  geom_errorbar(aes(x=Season, ymin = lci_66, ymax = uci_66), width = 0  ,  
                position = position_dodge(width = 0.8), size = 0.9 )+ 
  scale_color_viridis_d(option = "plasma" ) + 
  ylab(bquote('Seasonal AN'[Delta*degree*C==.(i)](per~100~k~2010~pop)))+ #paste0("Total Attributable Deaths\nRelative to Hindcast at ",i ," Δ°C" ))  + 
    theme_classic() +
  facet_wrap( vars(`Historical period`))

print(ci_season_rr_2010100k_tp_facet)
 
 
ggsave( paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, "_an_season_ci_2010100k_tp_facet.pdf" ), plot = ci_season_rr_2010100k_tp_facet )
ggsave( paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, "_an_season_ci_2010100k_tp_facet.png"), plot = ci_season_rr_2010100k_tp_facet ) }
## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

Last historical period based seasonal changes

Plot just the estimates from the last historical period to show changes by season

for (i in 2:6) {
  an_season_ci_2010100k_tp4 <- city_ci_sums_season %>% 
  filter(season %in% c(1:4) & grepl("an_total", tm_var)& time_period == 4 & dT_c == i) %>% # select out just the seasonal values and rr
  select( -value  ) %>% rename(value = an_pop2010_100k) %>% 
  pivot_wider()  %>% # select out just the seasonal values and rr %>%
  ggplot(aes(x = Season, y = mean_val, color = Cluster, group = Cluster)) +
  geom_point(   position = position_dodge(width = 0.8) ) + 
  geom_errorbar(aes(x=Season, ymin = lci_66, ymax = uci_66),
                width = 0  ,  position = position_dodge(width = 0.8) )+ 
  scale_color_viridis_d(option = "plasma" ) + 
  ylab( bquote('Seasonal AN'[Delta*degree*C==.(i)](per~100~k~2010~pop)))+ #paste0("Total Attributable Deaths Relative to Hindcast\nRelative to Hindcast at ",i ," Δ°C" )) +
    theme_classic() 

print(an_season_ci_2010100k_tp4)
 
 
ggsave(paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, 
              "_an_season_ci_2010100k_tp4.pdf" ), plot = an_season_ci_2010100k_tp4 )
ggsave(paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i,
              "_an_season_ci_2010100k_tp4.png"), plot = an_season_ci_2010100k_tp4 )  }
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

Plot the differences by historical period across clusters.

for (i in 2:6) {
  an_yr_ci_2010100k <- city_ci_sums_year %>% 
  filter(  grepl("an_total", tm_var) & dT_c == i) %>% # select out just the seasonal values and an
  select( -value  ) %>% rename(value = an_pop2010_100k) %>% 
  pivot_wider()  %>%  
  ggplot(aes(group =  `Historical period` , 
             x = Cluster, y = mean_val, color = `Historical period`)) +
  geom_point(position = position_dodge(width = 0.8) ) + 
  geom_errorbar(aes(x=Cluster, ymin = lci_66, ymax = uci_66), width = 0  , 
                position = position_dodge(width = 0.8) )+ 
  scale_color_viridis_d(  ) + 
  ylab( bquote('Yearly AN '[Delta*degree*C==.(i)](per~100~k~2010~pop)))  +# paste0("Yearly Temperature Attributable Deaths\nRelative to Hindcast at ",i ," Δ°C" )) 
  theme_classic() 
print(an_yr_ci_2010100k )
 
 
ggsave(paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, 
              "_an_yr_ci_2010100k.pdf" ), plot = an_yr_ci_2010100k )
ggsave(paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, 
              "_an_yr_ci_2010100k.png"), plot = an_yr_ci_2010100k )  }
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image
## Saving 7 x 5 in image

Plot the death attributable to mortality against the temperature ranges. The coldest are below the loweset exterior knot, the hottest above the highest exterior knot. The Cold are between the lowest exterior and the lowest interior. The Mid are between the lowest interior and the highest interior. The Hot are between the highest interiror and the highest exterior.

for (i in 2:6) {# i = 2
an_temp_ci_2010100k_facet_tp <- city_ci_sums_tmprng %>%   
  filter(  grepl("an_total", tm_var) & dT_c == i) %>% # select out just the seasonal values and an
  select( -value  ) %>% rename(value = an_pop2010_100k) %>%
  pivot_wider()  %>% 
  ggplot(aes(x = `Range C`, y = mean_val, group = Cluster, color = Cluster)) +
  geom_point(   position = position_dodge(width = 0.8) ) + 
  geom_errorbar(aes(x=`Range C`, ymin = lci_66, ymax = uci_66), width = 0  ,  position = position_dodge(width = 0.8) )+ 
  ylab(bquote('AN'[Delta*degree*C==.(i)](per~100~k~2010~pop))  ) + #paste0("Total Attributable Death Relative to Hindcast\nRelative to Hindcast at ",i ," Δ°C" )
  theme_classic() +
  scale_color_viridis_d(option = "plasma" ) + 
  facet_wrap( vars(`Historical period`))
expression("Temperature " ( degree~C))

print(an_temp_ci_2010100k_facet_tp )
 
ggsave(paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, 
              "_an_temp_ci_2010100k_facet_tp.pdf" ), plot = an_temp_ci_2010100k_facet_tp )
ggsave(paste0(here("output/y4-final-manuscript-output/2020-02-24-QA-Manu-fig-tab-PROJ/"), "/dT_C", i, 
              "_an_temp_ci_2010100k_facet_tp.png"), plot = an_temp_ci_2010100k_facet_tp )    }
## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

## Saving 7 x 4 in image
## Saving 7 x 4 in image

devtools::session_info()
## - Session info ---------------------------------------------------------------
##  setting  value                       
##  version  R version 3.6.2 (2019-12-12)
##  os       Windows 10 x64              
##  system   x86_64, mingw32             
##  ui       RTerm                       
##  language (EN)                        
##  collate  English_United States.1252  
##  ctype    English_United States.1252  
##  tz       America/Denver              
##  date     2020-02-26                  
## 
## - Packages -------------------------------------------------------------------
##  package     * version date       lib source        
##  assertthat    0.2.1   2019-03-21 [1] CRAN (R 3.6.0)
##  backports     1.1.5   2019-10-02 [1] CRAN (R 3.6.1)
##  broom         0.5.4   2020-01-27 [1] CRAN (R 3.6.2)
##  callr         3.4.2   2020-02-12 [1] CRAN (R 3.6.2)
##  cellranger    1.1.0   2016-07-27 [1] CRAN (R 3.6.0)
##  cli           2.0.1   2020-01-08 [1] CRAN (R 3.6.2)
##  colorspace    1.4-1   2019-03-18 [1] CRAN (R 3.6.0)
##  crayon        1.3.4   2017-09-16 [1] CRAN (R 3.6.0)
##  DBI           1.1.0   2019-12-15 [1] CRAN (R 3.6.2)
##  dbplyr        1.4.2   2019-06-17 [1] CRAN (R 3.6.2)
##  desc          1.2.0   2018-05-01 [1] CRAN (R 3.6.1)
##  devtools      2.2.2   2020-02-17 [1] CRAN (R 3.6.2)
##  digest        0.6.24  2020-02-12 [1] CRAN (R 3.6.2)
##  dlnm        * 2.3.9   2019-03-11 [1] CRAN (R 3.6.0)
##  dplyr       * 0.8.4   2020-01-31 [1] CRAN (R 3.6.2)
##  ellipsis      0.3.0   2019-09-20 [1] CRAN (R 3.6.2)
##  evaluate      0.14    2019-05-28 [1] CRAN (R 3.6.2)
##  fansi         0.4.1   2020-01-08 [1] CRAN (R 3.6.2)
##  farver        2.0.3   2020-01-16 [1] CRAN (R 3.6.2)
##  forcats     * 0.4.0   2019-02-17 [1] CRAN (R 3.6.0)
##  fs            1.3.1   2019-05-06 [1] CRAN (R 3.6.0)
##  generics      0.0.2   2018-11-29 [1] CRAN (R 3.6.0)
##  ggplot2     * 3.2.1   2019-08-10 [1] CRAN (R 3.6.2)
##  glue          1.3.1   2019-03-12 [1] CRAN (R 3.6.0)
##  gtable        0.3.0   2019-03-25 [1] CRAN (R 3.6.0)
##  haven         2.2.0   2019-11-08 [1] CRAN (R 3.6.2)
##  here        * 0.1     2017-05-28 [1] CRAN (R 3.6.0)
##  highr         0.8     2019-03-20 [1] CRAN (R 3.6.0)
##  hms           0.5.3   2020-01-08 [1] CRAN (R 3.6.2)
##  htmltools     0.4.0   2019-10-04 [1] CRAN (R 3.6.2)
##  httr          1.4.1   2019-08-05 [1] CRAN (R 3.6.2)
##  jsonlite      1.6.1   2020-02-02 [1] CRAN (R 3.6.2)
##  kableExtra  * 1.1.0   2019-03-16 [1] CRAN (R 3.6.2)
##  knitr       * 1.28    2020-02-06 [1] CRAN (R 3.6.2)
##  labeling      0.3     2014-08-23 [1] CRAN (R 3.6.0)
##  lattice       0.20-38 2018-11-04 [2] CRAN (R 3.6.2)
##  lazyeval      0.2.2   2019-03-15 [1] CRAN (R 3.6.0)
##  lifecycle     0.1.0   2019-08-01 [1] CRAN (R 3.6.2)
##  lubridate   * 1.7.4   2018-04-11 [1] CRAN (R 3.6.0)
##  magrittr      1.5     2014-11-22 [1] CRAN (R 3.6.0)
##  Matrix        1.2-18  2019-11-27 [2] CRAN (R 3.6.2)
##  memoise       1.1.0   2017-04-21 [1] CRAN (R 3.6.1)
##  mgcv          1.8-31  2019-11-09 [1] CRAN (R 3.6.2)
##  mixmeta       1.0.7   2019-12-09 [1] CRAN (R 3.6.2)
##  modelr        0.1.5   2019-08-08 [1] CRAN (R 3.6.2)
##  munsell       0.5.0   2018-06-12 [1] CRAN (R 3.6.0)
##  mvmeta      * 1.0.3   2019-12-10 [1] CRAN (R 3.6.2)
##  nlme          3.1-144 2020-02-06 [1] CRAN (R 3.6.2)
##  pillar        1.4.3   2019-12-20 [1] CRAN (R 3.6.2)
##  pkgbuild      1.0.6   2019-10-09 [1] CRAN (R 3.6.2)
##  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 3.6.2)
##  pkgload       1.0.2   2018-10-29 [1] CRAN (R 3.6.1)
##  prettyunits   1.1.1   2020-01-24 [1] CRAN (R 3.6.2)
##  processx      3.4.2   2020-02-09 [1] CRAN (R 3.6.2)
##  ps            1.3.2   2020-02-13 [1] CRAN (R 3.6.2)
##  purrr       * 0.3.3   2019-10-18 [1] CRAN (R 3.6.2)
##  R6            2.4.1   2019-11-12 [1] CRAN (R 3.6.2)
##  Rcpp          1.0.3   2019-11-08 [1] CRAN (R 3.6.2)
##  readr       * 1.3.1   2018-12-21 [1] CRAN (R 3.6.0)
##  readxl      * 1.3.1   2019-03-13 [1] CRAN (R 3.6.0)
##  remotes       2.1.1   2020-02-15 [1] CRAN (R 3.6.2)
##  reprex        0.3.0   2019-05-16 [1] CRAN (R 3.6.0)
##  rlang         0.4.4   2020-01-28 [1] CRAN (R 3.6.2)
##  rmarkdown     2.1     2020-01-20 [1] CRAN (R 3.6.2)
##  rprojroot     1.3-2   2018-01-03 [1] CRAN (R 3.6.0)
##  rstudioapi    0.11    2020-02-07 [1] CRAN (R 3.6.2)
##  rvest         0.3.5   2019-11-08 [1] CRAN (R 3.6.2)
##  scales        1.1.0   2019-11-18 [1] CRAN (R 3.6.2)
##  sessioninfo   1.1.1   2018-11-05 [1] CRAN (R 3.6.1)
##  stringi       1.4.3   2019-03-12 [1] CRAN (R 3.6.0)
##  stringr     * 1.4.0   2019-02-10 [1] CRAN (R 3.6.0)
##  testthat      2.3.1   2019-12-01 [1] CRAN (R 3.6.2)
##  tibble      * 2.1.3   2019-06-06 [1] CRAN (R 3.6.2)
##  tidyr       * 1.0.2   2020-01-24 [1] CRAN (R 3.6.2)
##  tidyselect    1.0.0   2020-01-27 [1] CRAN (R 3.6.2)
##  tidyverse   * 1.3.0   2019-11-21 [1] CRAN (R 3.6.2)
##  tsModel       0.6     2013-06-24 [1] CRAN (R 3.6.0)
##  usethis       1.5.1   2019-07-04 [1] CRAN (R 3.6.1)
##  vctrs         0.2.3   2020-02-20 [1] CRAN (R 3.6.2)
##  viridisLite   0.3.0   2018-02-01 [1] CRAN (R 3.6.0)
##  webshot       0.5.2   2019-11-22 [1] CRAN (R 3.6.2)
##  withr         2.1.2   2018-03-15 [1] CRAN (R 3.6.0)
##  xfun          0.12    2020-01-13 [1] CRAN (R 3.6.2)
##  xml2          1.2.2   2019-08-09 [1] CRAN (R 3.6.2)
##  yaml          2.2.1   2020-02-01 [1] CRAN (R 3.6.2)
## 
## [1] C:/Users/LayC/R.Lib
## [2] C:/Program Files/R/R-3.6.2/library