Important Note

Data Sources certificates_2020: https://epc.opendatacommunities.org/domestic/search?address=&postcode=&local-authority=&constituency=&uprn=&from-month=1&from-year=2020&to-month=12&to-year=2020 authority_codes: https://public.opendatasoft.com/explore/dataset/georef-united-kingdom-local-authority-district/export/?disjunctive.ctry_code&disjunctive.ctry_name&disjunctive.rgn_code&disjunctive.rgn_name&disjunctive.ctyua_code&disjunctive.ctyua_name&disjunctive.lad_code&disjunctive.lad_name incomes_2020: https://www.ons.gov.uk/employmentandlabourmarket/peopleinwork/earningsandworkinghours/datasets/smallareaincomeestimatesformiddlelayersuperoutputareasenglandandwales housing_age: https://www.datadaptive.com/docs/Age%20Profiles%20of%20Housing%20Stock%20-%20User%20Guide.pdf govt_spending_per_head: https://www.gov.uk/government/statistics/country-and-regional-analysis-2020/country-and-regional-analysis-november-2020 average_house_size: https://assets.publishing.service.gov.uk/media/5f047a01d3bf7f2be8350262/Size_of_English_Homes_Fact_Sheet_EHS_2018.pdf size: https://www.ons.gov.uk/datasets/TS050/editions/2021/versions/1 fiscal_spending: https://www.gov.uk/government/statistics/dluhc-enabled-spend-statistics-2021-22-and-2022-23

# Group by local authority (smallest level of governance) and calculate averages
ratings_county <- certificates_2020 %>% 
  group_by(LOCAL_AUTHORITY) %>% 
  summarise(
    avg_energyeff_current = mean(CURRENT_ENERGY_EFFICIENCY),
    avg_energyeff_potential = mean(POTENTIAL_ENERGY_EFFICIENCY),
    avg_envt_impact_current = mean(ENVIRONMENT_IMPACT_CURRENT),
    avg_envt_impact_potential = mean(ENVIRONMENT_IMPACT_POTENTIAL),
    avg_CO2_emissions_current = mean(CO2_EMISSIONS_CURRENT),
    avg_CO2_emissions_potential = mean(CO2_EMISSIONS_POTENTIAL),
    avg_lighting_cost_current = mean(LIGHTING_COST_CURRENT),
    avg_lighting_cost_potential = mean(LIGHTING_COST_POTENTIAL),
    avg_hot_water_cost_current = mean(HOT_WATER_COST_CURRENT),
    avg_hot_water_cost_potential = mean(HOT_WATER_COST_POTENTIAL),
    avg_floor_area = mean(TOTAL_FLOOR_AREA)
  ) %>%
  mutate(
    ENERGY_CONSUMPTION_DIFFERENTIAL = avg_energyeff_potential - avg_energyeff_current
  ) %>%
  arrange(desc(ENERGY_CONSUMPTION_DIFFERENTIAL))
# Merge data
merged_data <- merge(ratings_county, authority_codes, by.x = "LOCAL_AUTHORITY", by.y = "LOCAL_AUTHORITY") 
merged_data <- merge(merged_data, incomes_2020, by.x = "LOCAL_AUTHORITY", by.y = "Local authority code", all.x = TRUE)
merged_data <- merge(merged_data, housing_age, by.x = "LOCAL_AUTHORITY", by.y = "GEOCODE", all.x = TRUE)
# Plot thematic map- current energy efficiency
tmap_mode("view") # makes map interactive
your_data <- merged_data[!is.na(merged_data$Geo.Point), ]
your_data <- separate(your_data, Geo.Point, into = c("latitude_column", "longitude_column"), sep = ",")
coordinates <- st_as_sf(your_data, coords = c("longitude_column", "latitude_column"), crs = 4326)

breaks <- c(45,50,55,60,65,70)
colors <- c("blue", "cyan", "green","yellow","orange","red")

tm_shape(coordinates) +
  tm_bubbles(size = "avg_energyeff_current", 
             col = "avg_energyeff_current",
             breaks = breaks,
             palette = colors,
             title.size = 1, title.col = "black", title = "Avg Energy Efficiency Current") +
  tm_layout(legend.position = c("right", "bottom"))
# Plot thematic map- potential energy efficiency 
your_data <- merged_data[!is.na(merged_data$Geo.Point), ]
your_data <- separate(your_data, Geo.Point, into = c("latitude_column", "longitude_column"), sep = ",")
coordinates <- st_as_sf(your_data, coords = c("longitude_column", "latitude_column"), crs = 4326)

breaks <- c(70,75,80,85,90,95)
colors <- c("blue", "cyan", "green","yellow","orange","red")

tm_shape(coordinates) +
  tm_bubbles(size = "avg_energyeff_potential", 
             col = "avg_energyeff_potential",
             breaks = breaks,
             palette = colors,
             title.size = 1, title.col = "black", title = "Avg Energy Efficiency Potential") +
  tm_layout(legend.position = c("right", "bottom"))
# Plot thematic map-  energy efficiency differential
your_data <- merged_data[!is.na(merged_data$Geo.Point), ]
your_data <- separate(your_data, Geo.Point, into = c("latitude_column", "longitude_column"), sep = ",")
coordinates <- st_as_sf(your_data, coords = c("longitude_column", "latitude_column"), crs = 4326)

breaks <- c(10,13,16,19,22,25)
colors <- c("violet", "blue", "green","yellow","orange","red")

tm_shape(coordinates) +
  tm_bubbles(size = "ENERGY_CONSUMPTION_DIFFERENTIAL", 
             col = "ENERGY_CONSUMPTION_DIFFERENTIAL",
             breaks = breaks,
             palette = colors,
             title.size = 1, title.col = "black", title = "Avg Energy Efficiency Differential (unmet capacity for more energy-efficent housing)") +
  tm_layout(legend.position = c("right", "bottom"))
merged_data_region <- merged_data %>% 
  group_by(Region.Name) %>% 
  summarise(
    avg_energyeff_current_region = mean(avg_energyeff_current),
    avg_energyeff_potential_region = mean(avg_energyeff_potential),
    avg_envt_impact_current_region = mean(avg_envt_impact_current),
    avg_envt_impact_potential_region = mean(avg_envt_impact_potential),
    avg_netincomes = mean(Disposable.Income),
    reg_BPPRE1900P = mean(BPPRE1900P),
    reg_BP190018P = mean(BP190018P),
    reg_BP191929P = mean(BP191929P),
    reg_BP193039P = mean(BP193039P),
    reg_BP194554P = mean(BP194554P),
    reg_BP195564P = mean(BP195564P),
    reg_BP196572P = mean(BP196572P),
    reg_BP197382P = mean(BP197382P),
    reg_BP198392P = mean(BP198392P),
    reg_BP199399P = mean(BP199399P),
    reg_BP200008P = mean(BP200008P),
    reg_BP200911P = mean(BP200911P),
    reg_BP201214P = mean(BP201214P),
    reg_BP201517P = mean(BP201517P),
    reg_BP201820P = mean(BP201820P)
  ) %>%
  mutate(
    ENERGY_CONSUMPTION_DIFFERENTIAL_REGION = avg_energyeff_potential_region - avg_energyeff_current_region
  ) %>%
  arrange(desc(ENERGY_CONSUMPTION_DIFFERENTIAL_REGION))


# Line chart
df_region <- merged_data_region %>%
  select(Region.Name, matches("^reg_BP"))

dfl_region <- melt(df_region, id.vars = "Region.Name")
dfl_region$period <- str_extract(dfl_region$variable, "\\d+")
dfl_region$start <- as.numeric(str_sub(dfl_region$period, 1, 4))
dfl_region$end <- as.numeric(paste0(str_sub(dfl_region$period, 1, 2), str_sub(dfl_region$period, 5, 6)))
dfl_region$start[dfl_region$end == 19] <- 1890
dfl_region$end[dfl_region$end == 19] <- 1899
dfl_region$mid <- dfl_region$start + (dfl_region$end - dfl_region$start) / 2

ggplot(data = dfl_region, aes(x = mid, y = value, group = Region.Name)) +
  geom_line(aes(color = Region.Name)) +
  labs(x = "Year of Construction", y = "% of Housing in Region")

dfl_region <- dfl_region %>%
  mutate(age_2020 = 2020 - mid)

dfl_summary_region <- dfl_region %>%
  group_by(Region.Name) %>%
  summarise(total_houses_region = sum(value),
            weighted_age_sum_region = sum(age_2020 * value)) %>%
  mutate(average_age_region = weighted_age_sum_region / total_houses_region)
ggplot(data = merged_data_region, aes(x = avg_netincomes, y = avg_energyeff_current_region, label = Region.Name)) +
  geom_point(color = "green", size = 3) +
  geom_text(aes(label = Region.Name), nudge_x = 0.5, nudge_y = 0.5, size = 3) +
  labs(title = "Energy Efficiency vs Average Net Income",
       x = "Average Net Income (GBP)",
       y = "Average Energy Efficiency (%)") +
  theme_minimal()

# Regression analysis
merged_data_authority <- merged_data %>% 
  group_by(LOCAL_AUTHORITY) %>% 
  summarise(
    avg_energyeff_current_aut = mean(avg_energyeff_current),
    avg_energyeff_potential_aut = mean(avg_energyeff_potential),
    avg_envt_impact_current_aut = mean(avg_envt_impact_current),
    avg_envt_impact_potential_aut = mean(avg_envt_impact_potential),
    avg_netincomes = mean(Disposable.Income),
    avg_CO2_emissions_current_aut = mean(avg_CO2_emissions_current),
    avg_CO2_emissions_potential_aut = mean(avg_CO2_emissions_current),
    avg_lighting_cost_current_aut = mean(avg_lighting_cost_current),
    avg_lighting_cost_potential_aut = mean(avg_lighting_cost_potential),
    avg_hot_water_cost_current_aut = mean(avg_hot_water_cost_current),
    avg_hot_water_cost_potential_aut = mean(avg_hot_water_cost_potential),
    avg_floor_area_aut = mean(avg_floor_area),
    reg_BPPRE1900P = mean(BPPRE1900P),
    reg_BP190018P = mean(BP190018P),
    reg_BP191929P = mean(BP191929P),
    reg_BP193039P = mean(BP193039P),
    reg_BP194554P = mean(BP194554P),
    reg_BP195564P = mean(BP195564P),
    reg_BP196572P = mean(BP196572P),
    reg_BP197382P = mean(BP197382P),
    reg_BP198392P = mean(BP198392P),
    reg_BP199399P = mean(BP199399P),
    reg_BP200008P = mean(BP200008P),
    reg_BP200911P = mean(BP200911P),
    reg_BP201214P = mean(BP201214P),
    reg_BP201517P = mean(BP201517P),
    reg_BP201820P = mean(BP201820P)
  ) %>%
  mutate(
    ENERGY_CONSUMPTION_DIFFERENTIAL_AUT = avg_energyeff_potential_aut - avg_energyeff_current_aut
  ) %>%
  arrange(desc(ENERGY_CONSUMPTION_DIFFERENTIAL_AUT))

df <- merged_data_authority %>%
  select(LOCAL_AUTHORITY, matches("^reg_BP"))

dfl <- melt(df, id.vars = "LOCAL_AUTHORITY")
dfl$period <- str_extract(dfl$variable, "\\d+")
dfl$start <- as.numeric(str_sub(dfl$period, 1, 4))
dfl$end <- as.numeric(paste0(str_sub(dfl$period, 1, 2), str_sub(dfl$period, 5, 6)))
dfl$start[dfl$end == 19] <- 1890
dfl$end[dfl$end == 19] <- 1899
dfl$mid <- dfl$start + (dfl$end - dfl$start) / 2


dfl <- dfl %>%
  mutate(age_2020 = 2020 - mid)

dfl_summary <- dfl %>%
  group_by(LOCAL_AUTHORITY) %>%
  summarise(total_houses = sum(value),
            weighted_age_sum = sum(age_2020 * value)) %>%
  mutate(average_age = weighted_age_sum / total_houses)

all_data <- merge(merged_data_authority, dfl_summary, by.x = "LOCAL_AUTHORITY", by.y = "LOCAL_AUTHORITY", all.x = TRUE)

library(jtools)
reg1 <- lm(avg_energyeff_current_aut ~ average_age, data = all_data)
reg2 <- lm(avg_energyeff_current_aut ~ average_age + avg_floor_area_aut, data = all_data)
reg3 <- lm(avg_energyeff_current_aut ~ average_age + avg_floor_area_aut + avg_lighting_cost_current_aut, data = all_data)
reg4 <- lm(avg_energyeff_current_aut ~ average_age + avg_floor_area_aut + avg_lighting_cost_current_aut + avg_hot_water_cost_current_aut, data = all_data)

export_summs(reg1, reg2, reg3, reg4, model.names = c("Age", "Age & Size ","Age & Size & Lighting Cost", "Age & Size & Lighting Cost & Hot Water Cost"))
AgeAge & Size Age & Size & Lighting CostAge & Size & Lighting Cost & Hot Water Cost
(Intercept)74.24 ***90.17 ***100.41 ***103.42 ***
(1.09)   (1.75)   (1.95)   (1.25)   
average_age-0.13 ***-0.17 ***-0.16 ***-0.14 ***
(0.02)   (0.02)   (0.01)   (0.01)   
avg_floor_area_aut       -0.15 ***0.08 ** 0.09 ***
       (0.01)   (0.03)   (0.02)   
avg_lighting_cost_current_aut              -0.38 ***-0.27 ***
              (0.04)   (0.03)   
avg_hot_water_cost_current_aut                     -0.10 ***
                     (0.00)   
N329       329       329       329       
R20.15    0.37    0.50    0.80    
*** p < 0.001; ** p < 0.01; * p < 0.05.