si.Rmd

co2eq_10_10 = round((co2eq %>% filter(eu_q_rank == 10))$value/(co2eq %>% filter(eu_q_rank == 1))$value,digits = 1)

## total per decile

exp_bottom_decile = round((exp %>% filter(eu_q_rank == 1))$value, digits = 1)

exp_top_decile = round((exp %>% filter(eu_q_rank == 10))$value, digits = 1)

energy_bottom_decile = round((energy %>% filter(eu_q_rank == 1))$value, digits = 1)

energy_top_decile = round((energy %>% filter(eu_q_rank == 10))$value, digits = 1)

co2eq_bottom_decile = round((co2eq %>% filter(eu_q_rank == 1))$value, digits = 1)

co2eq_top_decile = round((co2eq %>% filter(eu_q_rank == 10))$value, digits = 1)


## per adult equivalent per decile

aeu = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005,
         indicator == "total_adult_eq") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value),
            eu_ntile_name = first(eu_ntile_name)) 

exp_pae = exp %>%
  rename(total_fd_me = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_fd_e = (total_fd_me/value)*1000000000000)

fd_pae_bottom_decile = round((exp_pae %>% filter(eu_q_rank == 1))$pae_fd_e, digits = 0)

fd_pae_top_decile = round((exp_pae %>% filter(eu_q_rank == 10))$pae_fd_e, digits = 0)

energy_pae = energy %>%
  rename(total_energy_use_tj = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_energy_use_gj = (total_energy_use_tj/value)*1000000000)

energy_pae_bottom_decile = round((energy_pae %>% filter(eu_q_rank == 1))$pae_energy_use_gj, digits = 1)

energy_pae_top_decile = round((energy_pae %>% filter(eu_q_rank == 10))$pae_energy_use_gj, digits = 1)

co2eq_pae = co2eq %>%
  rename(total_co2eq_kg = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_co2eq_t = (total_co2eq_kg/value)*1000000)

co2eq_pae_bottom_decile = round((co2eq_pae %>% filter(eu_q_rank == 1))$pae_co2eq_t, digits = 1)

co2eq_pae_top_decile = round((co2eq_pae %>% filter(eu_q_rank == 10))$pae_co2eq_t, digits = 1)


## intensities

mean_energy_intens = dat_country_summary_by_eu_ntile %>%
                 filter(year == 2005) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me))

mean_energy_intens_bottom_decile = round((mean_energy_intens %>% filter(eu_q_rank == 1))$pe_energy_use_mj, digits = 1)

mean_energy_intens_top_decile = round((mean_energy_intens %>% filter(eu_q_rank == 10))$pe_energy_use_mj, digits = 1)

mean_co2eq_of_energy_intens = dat_country_summary_by_eu_ntile %>%
  filter(year == 2005) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj) %>%
  group_by(eu_q_rank) %>%
  summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj))

mean_co2eq_of_energy_intens_bottom_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 1))$intensity_e_c, digits = 1)

mean_co2eq_of_energy_intens_top_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 10))$intensity_e_c, digits = 1)

```

In 2005, total household final demand was `r exp_total` trn€, the energy footprint `r energy_total` EJ, and the carbon footprint `r co2eq_total` MtCO2eq. We estimated the 10:10 ratios at `r exp_10_10` for expenditure, `r energy_10_10` for the energy footprint, and `r co2eq_10_10` for the carbon footprint. Total expenditure ranged from `r exp_bottom_decile` trn€ to `r exp_top_decile` trn€ (or `r fd_pae_bottom_decile`€ to `r fd_pae_top_decile`€ per adult equivalent) across bottom and top decile, the energy footprint from `r energy_bottom_decile` EJ to `r energy_top_decile` EJ (or `r energy_pae_bottom_decile` GJ/ae to `r energy_pae_top_decile` GJ/ae), and the GHG footprint from `r co2eq_bottom_decile` MtCO2eq to `r co2eq_top_decile` MtCO2eq (or `r co2eq_pae_bottom_decile` tCO2eq/ae to `r co2eq_pae_top_decile` tCO2eq/ae). The average energy intensity of consumption was `r mean_energy_intens_bottom_decile` MJ/€ in the bottom decile and `r mean_energy_intens_top_decile` MJ/€ in the top decile. The GHG intensity of energy use was `r mean_co2eq_of_energy_intens_bottom_decile` gCO2eq/TJ in the bottom decile, and `r mean_co2eq_of_energy_intens_top_decile` gCO2eq/TJ in the top decile.

```{r figureS1-test, fig.width=12, fig.height=8, out.width="98%", fig.cap="**Figure S1: Expenditure and resource footprints and intensities across European expenditure deciles in 2005. Total expenditures (a), energy footprint (b), and GHG footprint (c) per decile. Energy intensity as energy footprint per expenditure (d), GHG intensity as GHG footprint per expenditure (e), and GHG intensity as GHG footprint per energy footprint (f).**"}

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') +
  plot_layout(guides = 'collect')  & 
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        legend.position = 'bottom',
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))
a

ggsave(here("analysis", "figures", "figureS1.pdf"), device=cairo_pdf)

```

## 2010 using main method, EXIOBASE industry-by-industry

```{r ntiles-total-2010-ixi}

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_top = p1 + p2 + p3

```

```{r ntiles-intensity-violin-2010-ixi}

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (gCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_bottom = p2 + p1 + p3

```

```{r values-in-text-2010-ixi}

# values in text

## inequality
exp = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

exp_total = (exp %>% summarise(value = sum(value)))$value

exp_10_10 = round((exp %>% filter(eu_q_rank == 10))$value/(exp %>% filter(eu_q_rank == 1))$value,digits = 1)

energy = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

energy_total = (energy %>% summarise(value = sum(value)))$value

energy_10_10 = round((energy %>% filter(eu_q_rank == 10))$value/(energy %>% filter(eu_q_rank == 1))$value,digits = 1)

co2eq = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name))

co2eq_total = (co2eq %>% summarise(value = sum(value)))$value

co2eq_10_10 = round((co2eq %>% filter(eu_q_rank == 10))$value/(co2eq %>% filter(eu_q_rank == 1))$value,digits = 1)

## total per decile

exp_bottom_decile = round((exp %>% filter(eu_q_rank == 1))$value, digits = 1)

exp_top_decile = round((exp %>% filter(eu_q_rank == 10))$value, digits = 1)

energy_bottom_decile = round((energy %>% filter(eu_q_rank == 1))$value, digits = 1)

energy_top_decile = round((energy %>% filter(eu_q_rank == 10))$value, digits = 1)

co2eq_bottom_decile = round((co2eq %>% filter(eu_q_rank == 1))$value, digits = 1)

co2eq_top_decile = round((co2eq %>% filter(eu_q_rank == 10))$value, digits = 1)


## per adult equivalent per decile

aeu = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010,
         indicator == "total_adult_eq") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value),
            eu_ntile_name = first(eu_ntile_name)) 

exp_pae = exp %>%
  rename(total_fd_me = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_fd_e = (total_fd_me/value)*1000000000000)

fd_pae_bottom_decile = round((exp_pae %>% filter(eu_q_rank == 1))$pae_fd_e, digits = 0)

fd_pae_top_decile = round((exp_pae %>% filter(eu_q_rank == 10))$pae_fd_e, digits = 0)

energy_pae = energy %>%
  rename(total_energy_use_tj = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_energy_use_gj = (total_energy_use_tj/value)*1000000000)

energy_pae_bottom_decile = round((energy_pae %>% filter(eu_q_rank == 1))$pae_energy_use_gj, digits = 1)

energy_pae_top_decile = round((energy_pae %>% filter(eu_q_rank == 10))$pae_energy_use_gj, digits = 1)

co2eq_pae = co2eq %>%
  rename(total_co2eq_kg = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_co2eq_t = (total_co2eq_kg/value)*1000000)

co2eq_pae_bottom_decile = round((co2eq_pae %>% filter(eu_q_rank == 1))$pae_co2eq_t, digits = 1)

co2eq_pae_top_decile = round((co2eq_pae %>% filter(eu_q_rank == 10))$pae_co2eq_t, digits = 1)


## intensities

mean_energy_intens = dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me))

mean_energy_intens_bottom_decile = round((mean_energy_intens %>% filter(eu_q_rank == 1))$pe_energy_use_mj, digits = 1)

mean_energy_intens_top_decile = round((mean_energy_intens %>% filter(eu_q_rank == 10))$pe_energy_use_mj, digits = 1)

mean_co2eq_of_energy_intens = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj) %>%
  group_by(eu_q_rank) %>%
  summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj))

mean_co2eq_of_energy_intens_bottom_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 1))$intensity_e_c, digits = 1)

mean_co2eq_of_energy_intens_top_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 10))$intensity_e_c, digits = 1)

```

In 2010, total household final demand was `r exp_total` trn€, the energy footprint `r energy_total` EJ, and the carbon footprint `r co2eq_total` MtCO2eq. We estimated the 10:10 ratios at `r exp_10_10` for expenditure, `r energy_10_10` for the energy footprint, and `r co2eq_10_10` for the carbon footprint. Total expenditure ranged from `r exp_bottom_decile` trn€ to `r exp_top_decile` trn€ (or `r fd_pae_bottom_decile`€ to `r fd_pae_top_decile`€ per adult equivalent) across bottom and top decile, the energy footprint from `r energy_bottom_decile` EJ to `r energy_top_decile` EJ (or `r energy_pae_bottom_decile` GJ/ae to `r energy_pae_top_decile` GJ/ae), and the GHG footprint from `r co2eq_bottom_decile` MtCO2eq to `r co2eq_top_decile` MtCO2eq (or `r co2eq_pae_bottom_decile` tCO2eq/ae to `r co2eq_pae_top_decile` tCO2eq/ae). The average energy intensity of consumption was `r mean_energy_intens_bottom_decile` MJ/€ in the bottom decile and `r mean_energy_intens_top_decile` MJ/€ in the top decile. The GHG intensity of energy use was `r mean_co2eq_of_energy_intens_bottom_decile` gCO2eq/TJ in the bottom decile, and `r mean_co2eq_of_energy_intens_top_decile` gCO2eq/TJ in the top decile.

```{r figureS2-test, fig.width=12, fig.height=8, out.width="98%", fig.cap="**Figure S2: Expenditure and resource footprints and intensities across European expenditure deciles in 2010. Total expenditures (a), energy footprint (b), and GHG footprint (c) per decile. Energy intensity as energy footprint per expenditure (d), GHG intensity as GHG footprint per expenditure (e), and GHG intensity as GHG footprint per energy footprint (f).**"}

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') +
  plot_layout(guides = 'collect')  & 
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        legend.position = 'bottom',
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

a

ggsave(here("analysis", "figures", "figureS2.pdf"), device=cairo_pdf)

```

## 2010 using main method, EXIOBASE product-by-product

```{r load-data2, include=FALSE}
# load data wrangling functions
source(here("analysis", "R", "si", "wrangler_functions_pxp.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

```

```{r ntiles-total-2010-pxp}

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 


p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_top = p1 + p2 + p3

```

```{r ntiles-intensity-violin-2010-pxp}

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (gCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_bottom = p2 + p1 + p3

```

```{r values-in-text-2010-pxp}

# values in text

## inequality
exp = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

exp_total = (exp %>% summarise(value = sum(value)))$value

exp_10_10 = round((exp %>% filter(eu_q_rank == 10))$value/(exp %>% filter(eu_q_rank == 1))$value,digits = 1)

energy = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

energy_total = (energy %>% summarise(value = sum(value)))$value

energy_10_10 = round((energy %>% filter(eu_q_rank == 10))$value/(energy %>% filter(eu_q_rank == 1))$value,digits = 1)

co2eq = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name))

co2eq_total = (co2eq %>% summarise(value = sum(value)))$value

co2eq_10_10 = round((co2eq %>% filter(eu_q_rank == 10))$value/(co2eq %>% filter(eu_q_rank == 1))$value,digits = 1)

## total per decile

exp_bottom_decile = round((exp %>% filter(eu_q_rank == 1))$value, digits = 1)

exp_top_decile = round((exp %>% filter(eu_q_rank == 10))$value, digits = 1)

energy_bottom_decile = round((energy %>% filter(eu_q_rank == 1))$value, digits = 1)

energy_top_decile = round((energy %>% filter(eu_q_rank == 10))$value, digits = 1)

co2eq_bottom_decile = round((co2eq %>% filter(eu_q_rank == 1))$value, digits = 1)

co2eq_top_decile = round((co2eq %>% filter(eu_q_rank == 10))$value, digits = 1)


## per adult equivalent per decile

aeu = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010,
         indicator == "total_adult_eq") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value),
            eu_ntile_name = first(eu_ntile_name)) 

exp_pae = exp %>%
  rename(total_fd_me = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_fd_e = (total_fd_me/value)*1000000000000)

fd_pae_bottom_decile = round((exp_pae %>% filter(eu_q_rank == 1))$pae_fd_e, digits = 0)

fd_pae_top_decile = round((exp_pae %>% filter(eu_q_rank == 10))$pae_fd_e, digits = 0)

energy_pae = energy %>%
  rename(total_energy_use_tj = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_energy_use_gj = (total_energy_use_tj/value)*1000000000)

energy_pae_bottom_decile = round((energy_pae %>% filter(eu_q_rank == 1))$pae_energy_use_gj, digits = 1)

energy_pae_top_decile = round((energy_pae %>% filter(eu_q_rank == 10))$pae_energy_use_gj, digits = 1)

co2eq_pae = co2eq %>%
  rename(total_co2eq_kg = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_co2eq_t = (total_co2eq_kg/value)*1000000)

co2eq_pae_bottom_decile = round((co2eq_pae %>% filter(eu_q_rank == 1))$pae_co2eq_t, digits = 1)

co2eq_pae_top_decile = round((co2eq_pae %>% filter(eu_q_rank == 10))$pae_co2eq_t, digits = 1)


## intensities

mean_energy_intens = dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me))

mean_energy_intens_bottom_decile = round((mean_energy_intens %>% filter(eu_q_rank == 1))$pe_energy_use_mj, digits = 1)

mean_energy_intens_top_decile = round((mean_energy_intens %>% filter(eu_q_rank == 10))$pe_energy_use_mj, digits = 1)

mean_co2eq_of_energy_intens = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj) %>%
  group_by(eu_q_rank) %>%
  summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj))

mean_co2eq_of_energy_intens_bottom_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 1))$intensity_e_c, digits = 1)

mean_co2eq_of_energy_intens_top_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 10))$intensity_e_c, digits = 1)

```

In 2010 but using the product-by-product version of EXIOBASE, total household final demand was `r exp_total` trn€, the energy footprint `r energy_total` EJ, and the GHG footprint `r co2eq_total` MtCO2eq. We estimated the 10:10 ratios at `r exp_10_10` for expenditure, `r energy_10_10` for the energy footprint, and `r co2eq_10_10` for the carbon footprint. Total expenditure ranged from `r exp_bottom_decile` trn€ to `r exp_top_decile` trn€ (or `r fd_pae_bottom_decile`€ to `r fd_pae_top_decile`€ per adult equivalent) across bottom and top decile, the energy footprint from `r energy_bottom_decile` EJ to `r energy_top_decile` EJ (or `r energy_pae_bottom_decile` GJ/ae to `r energy_pae_top_decile` GJ/ae), and the GHG footprint from `r co2eq_bottom_decile` MtCO2eq to `r co2eq_top_decile` MtCO2eq (or `r co2eq_pae_bottom_decile` tCO2eq/ae to `r co2eq_pae_top_decile` tCO2eq/ae). The average energy intensity of consumption was `r mean_energy_intens_bottom_decile` MJ/€ in the bottom decile and `r mean_energy_intens_top_decile` MJ/€ in the top decile. The GHG intensity of energy use was `r mean_co2eq_of_energy_intens_bottom_decile` gCO2eq/TJ in the bottom decile, and `r mean_co2eq_of_energy_intens_top_decile` gCO2eq/TJ in the top decile.

```{r figureS3-test, fig.width=12, fig.height=8, out.width="98%", fig.cap="**Figure S3: Expenditure and resource footprints and intensities across European expenditure deciles in 2010 (product-by-product EXIOBASE version). Total expenditures (a), energy footprint (b), and GHG footprint (c) per decile. Energy intensity as energy footprint per expenditure (d), GHG intensity as GHG footprint per expenditure (e), and GHG intensity as GHG footprint per energy footprint (f).**"}

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') +
  plot_layout(guides = 'collect')  & 
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        legend.position = 'bottom',
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

a

ggsave(here("analysis", "figures", "figureS3.pdf"), device=cairo_pdf)

```

## 2015 using alternative method, EXIOBASE industry-by-industry

There is relatively good agreement between the alternative method footprints and the EXIOBASE footprints, except in Bulgaria, where the alternative method footprint was around 3 times larger. Including Bulgaria, the alternative method footprints were on average 20% larger than the EXIOBASE footprints. With Bulgaria removed they were 10% larger on average. Eastern European countries especially had larger alternative method footprints due to high intensities in electricity production and hot water supply especially, and then more expenditure in CP045 multiplied by these intensities than the expenditure in EXIOBASE. This is why the figure shows such high footprints in the bottom European deciles compared to the method keeping EXIOBASE footprints the same (ms results).

```{r load-data3, include=FALSE}
# load data wrangling functions
source(here("analysis", "R", "si", "wrangler_functions_method2_ixi.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "coicop")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "coicop")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

```

```{r ntiles-total-2015-method2}

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_top = p1 + p2 + p3

```

```{r ntiles-intensity-violin-2015-method2}

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (gCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_bottom = p2 + p1 + p3

```

```{r figureS4-test, fig.width=12, fig.height=8, out.width="98%", fig.cap="**Figure S4: Expenditure and resource footprints and intensities across European expenditure deciles in 2015 (EXIOBASE industry-by-industry version, but using alternative method). Total expenditures (a), energy footprint (b), and GHG footprint (c) per decile. Energy intensity as energy footprint per expenditure (d), GHG intensity as GHG footprint per expenditure (e), and GHG intensity as GHG footprint per energy footprint (f).**"}

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') +
  plot_layout(guides = 'collect')  & 
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        legend.position = 'bottom',
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

a

ggsave(here("analysis", "figures", "figureS4.pdf"), device=cairo_pdf)

```

# Units of analysis

```{r load-data4, include=FALSE}
# load data wrangling functions
source(here("analysis", "R", "wrangler_functions.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

```

```{r }

# ae vs. per capita

dat_adult_eq = read_csv(here("analysis", "data", "derived", "si", "adult_eq_per_household.csv")) %>%
  filter(year==2015)

dat_tmp = dat_country_summary_by_cquint_and_euntile %>%
  filter(year == 2015) %>%
  left_join(dat_adult_eq %>%
              select(iso2, iso3, quint, adult_e_p_hh), by=c("iso2", "quint")) 

europe28 = c("AUT",
             "BEL",
             "BGR",
             "CYP",
             "CZE",
             "DEU",
             "DNK",
             "EST",
             "GRC",
             "ESP",
             "FIN",
             "FRA",
             "HRV",
             "HUN",
             "IRL",
             "LTU",
             "LVA",
             "MLT",
             "NLD",
             "NOR",
             "POL",
             "PRT",
             "ROU",
             "SWE",
             "SVN",
             "SVK",
             "TUR",
             "GBR")

population = wbstats::wb(country = europe28, 
                         indicator = "SP.POP.TOTL", 
                         startdate = 1990, enddate = 2017) %>%
  select(iso3 = iso3c, population = value, year = date) %>%
  mutate(year = parse_number(year))

dat_people_per_adult_eq = dat_tmp %>%
  ungroup() %>% 
  select(year, iso2, iso3, quint, total_adult_eq) %>%
  group_by(year, iso2, iso3) %>%
  summarise(total_adult_eq = sum(total_adult_eq)) %>%
  left_join(population, by = c("year","iso3")) %>%
  mutate(pop_per_ae = population/total_adult_eq) %>%
  ungroup() %>%
  select(iso2, pop_per_ae)

dat_tmp2 = dat_tmp %>%
  left_join(dat_people_per_adult_eq, by = "iso2") %>%
  mutate(pop_estimate_per_quintile = total_adult_eq*pop_per_ae) %>%
  ungroup() %>%
  select(iso2,iso3,quint,eu_q_rank,
                              total_adult_eq,
                              pop_estimate_per_quintile) %>%
  rename(total_population = pop_estimate_per_quintile)

eu_ntile_name = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")

a_tmp = dat_tmp2 %>%
  group_by(eu_q_rank) %>%
  summarise(total_population = sum(total_population)) %>%
  cbind(eu_ntile_name)

a = dat_tmp2 %>%
  group_by(eu_q_rank) %>%
  summarise(total_adult_eq = sum(total_adult_eq),
            total_population = sum(total_population)) %>%
  cbind(eu_ntile_name) %>%
  pivot_longer(cols = c(-eu_q_rank,-eu_ntile_name), names_to = "indicator", values_to="value") %>%
  mutate(indicator = dplyr::recode(indicator,
                                   "total_adult_eq" = "Adult equivalents",
                                   "total_population" = "Total population")) %>%
  mutate(value = value/1000000)

```

```{r values-in-text-ae-vs-pop}

ae_vs_pop = dat_tmp2 %>%
  summarise(total_adult_eq = sum(total_adult_eq),
            total_population = sum(total_population)) 

ae_share_of_pop = ae_vs_pop$total_adult_eq/ae_vs_pop$total_population

```

Through the main paper we use household per adult equivalent as our unit of analysis, following the EUROSTAT HBS. This meant that we adjusted decarbonisation scenario final energy numbers from total per capita to household per adult equivalent to better compare them with our environmental footprint estimates. We adjusted them for 1) the household share of the total footprint, and 2) the adult equivalent share of the total population.

As a numerical example, we adjust a total final energy of 53 GJ per capita from the LED scenario (Grubler et al. (2018) [@grubler_low_2018]), first by the household share of the total European energy footprint in 2015 (around 0.62, calculated in EXIOBASE), and then the share of total adult equivalents in the total European population in 2015 (also around `r ae_share_of_pop`, calculated using the EUROSTAT HBS, number of households per country, and population data per country). A total final energy of 53 GJ/capita is therefore adjusted to a household final energy of 53 GJ/adult equivalent in Europe ((53 total GJ/capita * 0.62 household share of total footprint)/0.62 adult equivalent share of total population = 53 household GJ/adult equivalent).

Our European expenditure deciles were constructed having the exact same number of adult equivalents per decile. When comparing with external per capita numbers, however, there are not the same number of population per decile because of differences in non-adult-equivalent-normalized people per household between income quintiles per country, and between countries. In Figure Sx we show an estimate of population per European expenditure decile. We use this to re-estimate our energy footprint per European expenditure decile in per capita terms, and then re-create Figure 5 from the main paper (Figure Sxx below) in per capita terms. 

Because we have the number of adult equivalents per country, and the total population per country, we could use both to calculate the total population per adult equivalent ratio for each country. We applied these ratio to the adult equivalents from different countries making up each European expenditure decile, as so we could estimate total population per European expenditure decile taking into account differences in household per capita between countries, but not between income quintiles within each country. Figure Sx shows these population estimates per European expenditure decile. 

```{r figureSx, fig.cap="**Figure S5:**"}

ggplot(a, aes(x=factor(eu_ntile_name), y=value)) +
  geom_col() +
  ylab("number (in millions)") +
  xlab("") +
  facet_wrap(~indicator, scales="free_y") +
  theme_ipsum() +
  theme(axis.text.x = element_text(angle = 90)) 

ggsave(here("analysis", "figures", "figureSx.pdf"), device=cairo_pdf)

```

## Figure 5 from manuscript in household final energy per capita

Here we re-create Figure 5 from the main paper after estimating our energy footprint per European expenditure decile in per capita terms (using the population per decile estimates from above), instead of per adult equivalent terms. Now we only adjust the decarbonisation scenario final energy numbers from total GJ per capita to household GJ per capita, so 53 total GJ/capita becomes: 53 * 0.62 = 33 household GJ/capita. 

```{r }

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>% 
  mutate(intensity_energy = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank, 
         intensity_energy) %>%
  filter(eu_q_rank == 10) %>%
  select(-eu_q_rank)

pdat_final_demand = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  mutate(total_energy_use_tj_diff = total_energy_use_tj-total_energy_use_tj_new) %>%
  select(eu_q_rank,total_energy_use_tj_new) %>%
  group_by(eu_q_rank) %>%
  summarise(total_energy_use_tj_new = sum(total_energy_use_tj_new)) %>%
  left_join(a_tmp, by = "eu_q_rank") %>%
  mutate(pc_energy_use_tj = total_energy_use_tj_new/total_population,
         pc_energy_use_gj = pc_energy_use_tj*1000)

df_energy_deciles = pdat_final_demand %>%
  select(eu_q_rank, pc_energy_use_gj)

ineq_curr = df_energy_deciles$pc_energy_use_gj[10]/df_energy_deciles$pc_energy_use_gj[1]

```

```{r }

library(readxl)

df_scenario_info = read_excel(here("analysis/data/raw/scenarios_additions.xlsx"), sheet="overview") %>%
  select(scenario, fe_gj_pc = final_energy_gj_per_capita_2050,
         ccs_required = primary_energy_fossil_w_ccs2050_ej,
         description) %>%
  arrange(fe_gj_pc) %>%
  mutate(fe_gj_pc = fe_gj_pc*0.62,
         fe_gj_pc = round(fe_gj_pc),
         ccs_required = round(ccs_required)) 

mea = c(min(df_scenario_info$fe_gj_pc),max(df_scenario_info$fe_gj_pc))
mer = c(min(df_scenario_info$fe_gj_pc),53) #before was 33: multiplied 53 by 0.62 and rounded up

c_mean_mea = round((17+mea[2])/2) #multiplied 27 by 0.62 and rounded up
c_mean_mer = round((mer[1]+mer[2])/2)

```

```{r , eval = FALSE}

# run once to save file. If changing scenarios included or input data, need to re-run and save 'scenarios_fine.rds'

# vectorized function that returns scaled quantiles given