An Environmental Justice Analysis of Gas Leaks4.Rmd

---
title: "An Environmental Justice Analysis of Distribution-Level Natural Gas Leaks in Massachusetts, USA"
author:
date:
output: 
  bookdown::pdf_document2:
    latex_engine: xelatex
    keep_tex: yes
    toc: false
    number_sections: true
    fig_caption: yes        
    includes:  
      in_header: my_header.tex
link-citations: yes
bibliography: GasLeaksEJBIB.bib
csl: chicago-author-date.csl
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=FALSE, message=FALSE, warning=FALSE)
```
**Authors:**
Marcos Luna[^a][^*] and Dominic Nicholas[^b]

[^a]: Geography and Sustainability Department, Salem State University, 352 Lafayette Street, Salem, Massachusetts, USA 01970. mluna@salemstate.edu. 
[^b]: HEET (Home Energy Efficiency Team), 21 Acorn Street, Cambridge, Massachusetts, USA 02139. dominic.nicholas@heet.org.
[^*]: Corresponding author.


# Abstract {-}
A growing body of research shows that natural gas leaks at the distribution level are much more common and extensive than previously thought. Although scholars and advocates have raised alarms about the climate change and economic significance of these leaks, there has been little consideration of the problem from an environmental justice perspective. Using recently available high resolution leak data, this analysis of natural gas leaks across the state of Massachusetts shows that People of Color, limited English speaking households, renters, lower income residents, and adults with lower levels of education are disproportionately exposed to natural gas leaks and that their leaks take longer to repair, as compared to the general population, and particularly as compared to White residents and to homeowners. This pattern is evident for all leaks in the state, for leaks disaggregated by leak class or grade, and for leaks disaggregated by utility. This analysis shows that natural gas leaks are an environmental justice issue warranting further study and policy attention.

keywords: natural gas distribution, utilities, leaks, environmental justice, equity, race

\pagebreak
# Introduction

From the wellhead to the gas appliances in homes, researchers increasingly find that natural gas leaks are more common and extensive than previously understood or reported. These leaks have outsized climate impacts because methane, the primary constituent of natural gas, has more than 80 times the global warming potential of carbon dioxide over the first 20 years that it reaches the atmosphere. Lost methane may account for one quarter of current warming [@EDF2021]. Locally, gas leaks pose explosive fire hazards, contribute to degradation of indoor and outdoor air quality, and kill street trees. Moreover, the loss of natural gas represents an economic burden both to consumers and to distributors. Despite rapid growth in studies about the extent of gas leaks and methods to detect them, research on gas leaks from an equity perspective has been rare and inconclusive [@PHILLIPS20131; @Scott_Scott_2019]. As with many other forms of environmental and energy challenges, there is reason to suspect that the experience of natural gas leaks is neither equal nor equitable. 

This paper assesses the degree to which  gas leaks on distribution-level mains and services to buildings are distributed inequitably and how  repair of these leaks varies between communities across the state of Massachusetts, USA. To perform this assessment, we utilize newly available street-level records of gas leaks reported by natural gas utilities which have been geocoded by the Home Energy Efficiency Team (HEET), a non-profit environmental advocacy organization based in Cambridge, Massachusetts. The gas leaks records reveal not only detailed information about individual leak locations, but also their first reporting date, whether or when they were repaired, and their leak class or grade (i.e. hazardousness). In this paper, leak class and grade are used synonymously. We compare the relative frequency of these leaks across communities at varying geographic scales using American Community Survey census data to determine whether there are systematic differences in exposure or timeliness of repair for different geographic communities and population subgroups. 

This high resolution, statewide assessment of gas leaks exposure and response applies a much-needed environmental justice framework to the problem of natural gas leaks. It advances prior work on gas leaks occurrence which has been limited to independent analyses of specific municipalities [@PHILLIPS20131], or which has been conducted at coarse geographic scales due to the lack of detailed data below regional or utility-scale reporting [@Scott_Scott_2019]. This analysis draws on a robust body of environmental justice research and practice to demonstrate the application of an equity analysis to the problem of distribution-level gas leaks. The results validate concerns that exposure and response to natural gas leaks are inequitable. These findings affirm the value to regulators, utilities, and communities of adopting an environmental justice perspective to reveal and understand how inequities are created and perpetuated in natural gas distribution. This study also highlights the critical importance of detailed and transparent reporting by utilities and their regulators to support environmental justice. 


# Background
## Gas service provision in Massachusetts

Natural gas consumption in Massachusetts has grown dramatically in the last three decades – over 200% since the early 1980s – displacing coal, oil, and nuclear power for electricity generation, and wood and fuel oil for home heating [@EIA_2021]. The growth in natural gas consumption over this time period is variously attributed to the higher price and volatility of fuel oil for home heating [@Joskow2013], changes in federal regulation of the natural gas industry [@APGA2021; @NERC2011], the dramatic expansion of access to domestic shale gas through hydrologic fracking [@Logan2012; @NERC2011], and state-level regulation encouraging electric utilities to switch to cleaner and less carbon-intensive fuel sources (particularly compared to oil and coal) [@C2ES].At the end of 2019, approximately two-thirds of electricity generation in the state was fueled by natural gas, with the remainder made up primarily by hydroelectricity, biomass, and various sources of renewable energy [@EIA_2020].  However, households are the single largest consumers of natural gas, representing over 30% of total consumption for the state, most of which is used for hot water and space heating [@EIA_2015; @EIA_2019]. More than half (52%) of all residences in Massachusetts (1.4 million occupied housing units) rely on utility-supplied, piped natural gas as their primary source of home heating, with the remainder reliant primarily on fuel oil (24%) or electricity (17%) [@Census_2019].
Despite aggressive marketing by utilities and the natural gas trade association about the economic and environmental benefits of natural gas [@AGA2021; @Hall_2021; @Leber_2021], its appeal has been undermined by revelations about the extent of gas leaks, high profile disasters, and initiatives by the state and municipalities to dramatically and rapidly reduce greenhouse gas emissions [@Borunda2020; @Roberts_2019; @Volcovici2021]. 


## Crises of gas leakage

Leakage of natural gas from production through transportation and distribution has long been recognized, but a growing body of research has revealed that the extent of these losses is much greater than reported by the natural gas industry or estimated by regulators. Until 2009, the Environmental Protection Agency (EPA) relied on a joint EPA-industry assessment from 1996 that estimated about 1.4% of the methane in the natural gas system  was lost to the atmosphere, with incremental revisions of upstream emissions estimates in subsequent years [@Howarth_2014]. However, researchers using independent systems of analyses, including life cycle assessment and atmospheric measurements, found that methane emissions from the natural gas system were likely to be at least two to three times greater than government and industry estimates [@Brandt_Heath_2014; @Howarth2011; @McKain2015; @Plant_Kort_2019; @Inman_Grubert_2020]. Until recently, however, there was very little research or reporting allowing for disaggregation of lost gas or methane emissions at the distribution level [@Weller_Hamburg_2020]. 

In Massachusetts, @PHILLIPS20131 pioneered a portable, street-level method of gas leak detection using a vehicle-mounted cavity ring-down spectrometer. The results and the maps they released beginning in 2012 showed that Boston was suffused by excessive levels of methane along the 785 miles of road they surveyed [@Daley_2012]. They argued that these concentrations could not be explained by natural sources and must be from leaks in gas mains buried beneath the streets. Their methods were replicated in municipalities across the country, from Washington D.C. [@Jackson_Down_2014] and New York City [@Gallagher_Down_2015] to Los Angeles [@Hopkins_Kort_2016]. These studies confirmed that gas leaks are a ubiquitous problem in the distribution system. Despite the growth in such studies, however, most have continued to focus on methodologies for the detection and quantification of methane from natural gas distribution, with virtually no attention to its environmental health or social implications [@vonFischer_2017; @Cho_Ulrich_2020; @Keyes_Ridge_2020; @Weller_2018]. 

With over 6,000 miles of aging, leak-prone infrastructure [@DPU_2019], Massachusetts bears a disproportionate share of the country’s leak-prone gas pipelines [@Herdes_Campbell_2020]. Natural gas distribution systems are comprised of mains and services. Mains generally distribute gas into an area and services  deliver gas from the mains to homes or businesses. At the end of 2019, Massachusetts gas utilities reported over 4,600 miles of leak-prone cast iron or unprotected steel mains comprising 22% of their total mains mileage, compared to 6% for the USA. Similarly, over 171,000 services were made of the same leak-prone materials comprising over 13% of service lines in the state, compared to 4% for the USA [@PHMSA2019]. As one of the longest settled parts of the country, the infrastructure in Massachusetts reflects its age.  Approximately 14% of mains and 12% of services across Massachusetts were installed prior to 1960 with 14% and 8% respectively installed before 1940 [@PHMSA2019]. Beginning in the early 1990s, the Pipeline Hazards and Materials Safety Administration (PHMSA) began issuing alert bulletins on the need to replace cast iron pipe, and later unprotected steel pipe as well as any pipelines over 60 years of age, because of the higher rate of failure and leaks based on nationwide studies [@Herdes_Campbell_2020; @Vetter_Kuebel_2019]. In 2009, the Massachusetts Department of Public Utilities (DPU) started approving plans for accelerated replacement of leak-prone infrastructure [@DPU_2019]. However, the pace of repair was frequently criticized [@Norton2011]. In 2013, Congress released a study commissioned by Massachusetts Senator Edward Markey showing that natural gas customers across the country had paid \$20 billion over a decade (2000 - 2011) for gas that they had never received; over $1.5 billion in Massachusetts alone [@House_2013]. Massachusetts subsequently passed legislation requiring utilities to submit annual plans to repair or replace leak-prone natural gas infrastructure, adopt a consistent system of leak hazard classification with timelines for monitoring and repair, and crucially (for this analysis), to submit publicly-accessible reports of the status and locations of leaks to the DPU [@GasLeakAct]. This statewide program for public reporting of distribution-level gas leaks is the first such program in the USA that requires street-level geographic detail. However, utilities are not required to directly measure or report leak volume. 

The first gas leaks reports were filed with the DPU in 2015. HEET cleaned up and translated the data from these reports into accessible, online interactive maps showing the location, grade, and status of the leaks for hundreds of communities across the state. They released these maps to the public beginning in summer 2015 [@Chakrabarti_2015]. HEET has continued this service annually with each release of the DPU data, performing quality control, and becoming the de facto source for publicly accessible maps of utility-reported leaks in Massachusetts. The sustained visibility of gas leaks has become an invaluable resource for advocates to hold utilities, the DPU, and policy makers accountable and to address these leaks [@HEET_2017]. 


## Another gas disaster reveals systemic problems

In September 2018, concerns about the dangers of natural gas rose to the forefront again after a series of natural gas explosions and fires hit three communities in the Merrimack Valley of northeastern Massachusetts. As a result of overpressurization during gas line replacement by Columbia Gas, five homes were destroyed, and another 131 damaged by fire. People were injured and approximately 50,000 were forced to evacuate the area.  The City of Lawrence, a largely Hispanic, immigrant, and working-class community, was hardest hit by the crisis, accounting for all of the direct injuries and most of the property damage [@NTSB_2019; @Walters_2019; @Valencia_2020].

While the disaster was regarded as an unusual event, other investigations have highlighted systemic problems in how natural gas utilities across the state prioritize and handle gas leaks. Research released in 2016 showed that just 7% of leaks in the metro Boston distribution system accounted for 50% of gas leak emissions by volume [@Hendrick_Ackley_2016]. Moreover, these “super emitters” were just as likely to be Grade 3 leaks, which are not considered an explosion hazard and thus are not normally prioritized for repair. Massachusetts subsequently passed legislation requiring utilities to identify and prioritize Grade 3 leaks of “significant environmental impact” or SEIs [@EnergyDiversityAct]. Leak grades are described in more detail in section 3.1 below. Since 2017, HEET has collaborated with the state's largest gas utilities to help them adopt a reliable and proven new method for accurately identifying SEIs [@HEET2021]. However, HEET and others have found that utility repairs appear to eliminate all gas from the entire leak footprint only 30% of the time [@HEET2021; @Edwards2021]. 
An outside auditor contracted by the state in response to the Columbia Gas disaster found that most of the gas utilities in the state were moving too slowly in pipeline replacement to meet statutory obligations, were plagued by poor record keeping, inconsistent methods of leak monitoring and repair, and curiously, were improperly prioritizing repairs in suburban communities, rather than in more densely populated urban communities where more people were exposed and the risks higher [@Herdes_Campbell_2020].

Sustained attention to the problem of leaking natural gas for the last decade in Massachusetts has spurred awareness about the magnitude and extent of the problem. Most of this attention has been premised on concerns about greenhouse gas mitigation, cost recovery, and public safety. With very few exceptions [@Mascoop_2018; @Scott_Scott_2019], these conversations have not considered the distributional or procedural equity dimensions of the problem i.e. an environmental justice perspective.


## An Environmental Justice framework

Environmental justice is the principle that everyone has a right to a healthy and safe environment and to be treated fairly, with meaningful opportunity for participation in processes or decisions that affect them.  In practice, environmental justice is rooted in the recognition that specific communities and population groups – particularly People of Color and lower income people in the USA context – have been denied these rights [@Agyeman2016; @Baptista2019; @Bonorris2010; @Lado2017; @Lado2019; @Salazar_Clauson_2019; @USCC2016].

Since the 1980s, government reports, academic scholarship, and investigative journalism have repeatedly found that People of Color and low income communities across the country are disproportionately exposed or vulnerable to a wide range of environmental burdens. These communities are also subjected to procedural inequities through unequal treatment, enforcement, and opportunities for meaningful participation. Environmental justice research has revealed unequal exposure to hazardous waste [@Salazar_Clauson_2019], air pollution [@Colmer_Hardman_2020], water pollution [@Schaider_2019], noise pollution [@Collins_Nadybal_2020], traffic [@Pinto_Reichmuth_2019], and more recently, the climate change enhanced risks of flooding [@Hardy_Milligan_2017], excessive heat [@Hoffman_Shandas_2020], and extreme weather [@García-López_2018]. Conversely, these same communities and groups are disproportionately denied environmental benefits or amenities, such as access to greenspace [@Rigolon_Browning_2018], urban tree canopy coverage [@Schwarz_Fragkias_2015], lead-free housing [@Whitehead_LaToria_2019], and access to safe drinking water [@Fedinick_Taylor_2019]. Environmental justice scholarship has also shown that these communities suffer procedural inequities including unequal environmental compliance by private businesses [@McDonald_Jones_2018], unequal enforcement by government agencies [@Konisky_Reenock_2018, and less government aid for disaster recovery [@Howell2018]. 
These burdens do not occur singly or in isolation, but are overlapping, cumulative, and synergistic in their effects on health and wellbeing and even economic opportunity. 

The focus of an environmental justice perspective is therefore on: 

* those communities that have been, and continue to be, disproportionately burdened by pollution or other environmental insults, 
* those who are especially vulnerable to risks and threats, 
* those who have been treated differently in the enforcement of law or administration of public services,
* those who have been unfairly excluded from enjoying environmental benefits, and 
* those who have been denied a voice in decision making about their environments – “where they live, work, and play” (Robert Bullard quoted in Schweizer 1999).

We apply the environmental justice framework to gas leaks because the gas distribution infrastructure is a significant component of the urban environment in which people “live, work, and play.” Gas leaks have the potential to be a safety hazard to nearby residents, and the timely repair of all leaks represent important expressions of prioritization by utilities and government regulators for these urban environments and the people who live there. Systematic differences in the geographic distribution of leaks or their timely repair along the axes of race/ethnicity, class, or other dimensions of social stratification may be evidence of social bias or systemic inequity – an environmental injustice. Although this framework might overlap conceptually with scholarship on energy justice [@Jenkins2016] or even climate justice [@Simmons2020], we use the environmental justice framework because that is the language used by both community advocates and policy makers in Massachusetts and the USA in addressing a broad range of distributional and procedural environmental inequities, including those concerning energy infrastructure [@ClimateRoad; @WHEJAC2021; @ExecutiveOrder1994; @FERC_2021; @LopezNickerson_2021; @DOE_EJ2021; @Goldin2021; @Madaro2021; @Morales2021]. 

At its core, environmental justice is about the spatial distribution of environmental benefits and burdens [@Lee_2021]. This is because environmental phenomena, including energy infrastructure and services,  are inherently geographic. *Where* something happens or where resources and attention are directed is often closely related to *who* lives there. This analysis applies an environmental justice analysis to the phenomenon of natural gas leaks in Massachusetts using a geospatial approach to quantitatively describe and compare the experiences of different population groups with respect to both exposure to gas leaks and to the timeliness of their repair. 


# Data and methods
## Data

Gas leaks data for this analysis were acquired from HEET [@HEET2021a]. This data is derived from Annual Service Quality Reports (ASQR) that are submitted by investor-owned utilities to the DPU each February, which reflect gas leak activity for the previous calendar year [@HEET_2021; @GasLeakAct]. HEET geocoded the leaks based on their reported addresses or cross-street locations and performed extensive quality control on the results. HEET’s gas leak data contain all of the original data fields reported in the AQSRs, including addresses, date first reported, and leak class or grade of individual repaired and unrepaired gas leaks for each utility for one year – January 1 through December 31, 2019. A small number (86 out of 26,541) of the originally reported data were not usable in this analysis because of erroneous entries in the original reports. 

Gas leak data was disaggregated by repaired or unrepaired status and also by leak class or grade. Under Massachusetts regulations (220 CMR 114.04), leaks are classified by utilities according to the following classification system: 

* Grade 1 represents “an existing or probable hazard to persons or property.”  Requires “immediate commencement of repair and continuous action until the conditions are no longer hazardous, the source of the leak is eliminated, and permanent repairs have been completed.”
* Grade 2 is recognized as “nonhazardous to persons or property at the time of detection, but justifies scheduled repair based on probable future hazard.”  Must be repaired within 12 months from the date the leak was classified.  Must be reevaluated at least once every six months until eliminated.
* Grade 3 is recognized as “nonhazardous to persons or property at the time of detection and can be reasonably expected to remain nonhazardous.”  Must be reevaluated during the next scheduled survey, or within 12 months from the date last evaluated, whichever occurs first. Leaks classified on or after January 1, 2018 must be repaired or eliminated within eight years. 

As of 2018, gas utilities are also required to repair Grade 3 leaks of significant environmental impact (SEIs)  on an accelerated schedule, within one to three years of identification.  SEIs are not evaluated separately here due to the current variability in classification accuracy. Classification accuracy is expected to improve over time and should be reevaluated in the future as a separate category. 

In their ASQRs, all utilities reported leaks that were in a "repaired" or "unrepaired" status at the end of 2019.  Some of the utilities reported that some leaks were in an "eliminated" status, either due to pipeline replacement Gas System Enhancement Programs (GSEPs), or other reasons (for example, resurveys not finding gas at leak locations). These eliminated leaks were not consistently reported by all utilities in their ASQRs and may have been prioritized and repaired differently to the other "repaired" leaks. We chose not to include them in this analysis and recommend further research be focused on GSEPs and prioritization.

Population data was derived from the American Community Survey 5-year estimates for 2015-2019 at the Census Block Group level, Census Tract level, and the municipality level (i.e., Census County Subdivisions) via the tidycensus package version 0.11.4 in R (Walker and Herman 2021).   


## Methodology
### Population-weighted mean exposure

Exposure to gas leaks was measured as the population-weighted mean gas leaks density within the geographic unit of analysis (i.e., Block Group, Tract, or municipality). For example, gas leak density is calculated as the number of reported gas leaks in a Block Group divided by the area of the Block Group in square kilometers. Population-weighted mean gas leaks density is calculated by multiplying the gas leaks density of each Block Group and the respective population groups of the same Block Groups to get weighted values, from which a weighted mean is calculated. Only Block Groups falling within the natural gas utility service areas providing leak data were included in the analysis. 

A population-weighted mean is calculated according to the following formula:
$$ \overline{x}_w = \frac{\sum_{i=1}^{n} (x_i \times w_i)}{\sum_{i=1}^{n}w_i} $$
where $\overline{x}_w$ is the population-weighted mean of some variable, in this case gas leak density;

where ${x_i}$ represents an individual observation, in this case the leak density of a Block Group;

where ${w_i}$ is the weight associated with an observation, in this case the population of a block group.

Population-weighted means are a common approach in analyses of the distribution of environmental burdens and amenities, particularly when comparing between population subgroups [@Bell_Ebisu_2012; @Mikati_Benson_2018; @Richmond-Bryant_Mikati_2020]. The Environmental Protection Agency (EPA) uses population-weighted means as the primary way of expressing exposure differentials to a variety of environmental burdens in its nationwide EJSCREEN application [@EPA_2020]. One of the advantages of a population-weighted approach is that population estimates for Census statistical units (i.e., Block Groups or Tracts) are not correlated with population density. This is because Census units are constructed to have similar population sizes. In Massachusetts, over 90% of Block Groups have populations between 565 and 2,700 for the 2015-2019 period, which is comparable to the population range of Block Groups across the country. This means population weighting does not emphasize urban or high-density locations. If population-weighted gas leak density is higher in denser, more urban areas, this difference is not due to population weighting.  Instead, this difference likely reflects the underlying phenomenon itself or other drivers of the phenomenon [@EPA_2019]. 


### Relative exposure

To better compare differences between groups, results are presented as relative exposures. Relative exposure is calculated similarly to relative risk or risk ratio:
$$ RE = \frac{\overline{x}_w subgroup} {\overline{x}_w total pop} $$
where $RE$ is the relative exposure expressed as a ratio;

where $\overline{x}_w subgroup$ is the population-weighted mean exposure of a subgroup, such as People of Color;

where $\overline{x}_w total pop$ is the population-weighted mean exposure of the total or comparison population.

Relative exposure values greater than one indicate that the subgroup is more exposed than the general population. Values less than one indicate that the subgroup is less exposed than the general population. For example, a relative exposure of 1.5 would indicate that a group has an exposure that is 1.5 times that of, or 50% greater than, the total population. Risk or relative exposure ratios are common in both public health and in studies of environmental justice because of their simplicity in making comparisons and communicating differences in exposure between groups [@Debbage_2019; @Mikati_Benson_2018; @Harner2002]. 

For this analysis, relative exposures are calculated relative to the general population ACS estimates. For relative exposure by race, ethnicity, and income, total population estimates are used as the denominator. Total households are used as the denominator for household-level relative exposure. Total occupied housing units are used as the denominator for housing-unit level relative exposure, such as housing tenure or housing burden.


### Sensitivity Analyses

To assess potential bias from analytic choices, we performed sensitivity analyses to account for the impact of differences in definitions of vulnerable populations, uncertainty in ACS population estimates, and choices in scale or the geographic unit of analysis. 

This analysis looks at variations in exposure for a variety of demographic groups that environmental justice policy and research have identified as being especially vulnerable to environmental burdens, or deprived of environmental benefits, as a consequence of social or economic disadvantage, physical vulnerability, or historic and persistent discrimination and inequality. These groups include:

* People of Color (i.e., persons who are of Hispanic ethnicity or racially not White, including indigenous populations)
* In addition to the aggregate category of “People of Color”, non-Hispanic Asians, non-Hispanic Blacks, Hispanics, and non-Hispanic Whites were analyzed separately as well
* Low income persons (i.e., income less than twice the poverty line)
* Limited-English speaking households (i.e., households where no one over the age of 14 speaks English “very well”)
* Adults 25 years or older without a high school diploma
* Children under the age of 5
* Adults over the age of 64
* Disabled adults
* Renters
* Housing burdened households (i.e., households paying more than 30% of household income to rent or mortgage) 
* People living in Environmental Justice communities as defined by the 2021 *Climate Roadmap* legislation, which capture Census Block Groups with high proportions of People of Color, limited English speaking households, or lower income persons. Massachusetts Environmental Justice communities are represented by the prefix ‘MA’. 

Population numbers for this analysis were derived from the ACS 5-year estimates for the period 2015-2019. After the 2000 census, detailed data on demographic characteristics (e.g., income, education, housing) are only available from the ACS. While the decennial Census is based on a total enumeration of the population at one point in time, ACS estimates are based on a rolling sample of responses – about 3.5 million annually across the country – which are pooled from surveys compiled on a monthly basis across the year [@ACS_2020]. For small areas with populations less than 65,000, estimates are only available as pooled 5-year estimates. Because the ACS is based on a sample, estimates are subject to uncertainty or sampling error. This uncertainty is communicated through a Margin of Error (MOE) which accompanies each ACS estimate. For some population subgroups, such as racial and ethnic minorities, children, and renters, Census response rates can be low and uncertainty high [@OHare_2019]. Moreover, smaller geographic units, such as Block Groups, will also tend to have higher MOEs because they are derived from small samples. One approach to dealing with the high MOE of small area estimates is to aggregate groups (e.g., combining individual racial categories into one ‘People of Color’ category) or to use larger geographic units (e.g., Census Tracts rather than Census Block Groups). Both techniques reduce uncertainty for estimates, although this reduced uncertainty comes at the expense of demographic or spatial resolution. Both methods – aggregating groups and using larger geographic units - are applied here to consider the impacts of uncertainty and aggregation. 

A separate but related issue to the uncertainty of sampling and estimation is the Modifiable Areal Unit Problem (MAUP). MAUP refers to the fact that analytical results can be influenced by the geographic size or shape of units of analysis. In other words, different geographic definitions of a study area may result in different outcomes. The MAUP is a long-recognized challenge in geographic analysis [@Openshaw_Taylor_1979], but analysts are often constrained by data availability (e.g., administrative databases, sensor resolution, reporting conventions), as well as by lack of clear guidance or theory on an appropriate geographic unit of analysis [@Dark_Bram_2007]. Varying both the areal size and shape of the units of analysis can show whether results are dependent on choices for the scale and shape of the unit of analysis. Although cumbersome, this is recommended practice in environmental justice analyses [@Cutter_Holm_1996; @Baden_Noonan_Turaga_2007; @Karner2013]. To explore the impact of geographic scale and choice of unit of analysis, we replicated our analyses at the Census Block Group level, Census Tract level, and municipality level (i.e., Census County Subdivision) (see Figure \@ref(fig:MAUP) and Table \@ref(tab:MAUPtab)).

```{r MAUP, fig.align = "center", fig.cap="Comparison of scales and units of analysis: municipality (i.e., Census County Subdivision), Census Tract, and Census Block Group", echo=FALSE, message=FALSE}
# create a graphic to illustrate relationship of block groups, tracts, and county subdivisions
# library(tidyverse)
# library(tigris)
# library(tmap)
# library(sf)
# library(tmaptools)
# 
# ma_counties <- counties(state = "MA", cb = TRUE) %>%
#   filter(NAME == "Suffolk")
# # ma_cosub <- county_subdivisions(state = "MA", cb = TRUE, county = "025") %>%
# #   filter(NAME == "Chelsea") # NOT WORKING ALL OF SUDDEN ON 5/13/2021
# # Download file directly from https://www2.census.gov/geo/tiger/TIGER2019/
# ma_cosub <- st_read(dsn = "Data/tl_2019_25_cousub", layer = "tl_2019_25_cousub")%>%
#   filter(NAME == "Chelsea")
# ma_tracts <- tracts(state = "MA", cb = TRUE, county = "025")
# ma_blkgrps <- block_groups(state = "MA", cb = TRUE, county = "025")
# 
# 
# # isolate tracts in Chelsea
# ma_tracts_chelsea <- ma_tracts %>%
#   select(GEOID) %>%
#   st_centroid() %>%
#   filter(st_intersects(., ma_cosub, sparse = FALSE)) %>%
#   as.data.frame() %>%
#   right_join(ma_tracts, ., by = "GEOID")
# 
# # tm_shape(ma_tracts_chelsea) + tm_borders() +
# #   tm_shape(ma_tracts_chelsea) + tm_text(text = "NAME") #25025160101
# 
# # isolate tract
# ma_tracts_1 <- ma_tracts_chelsea %>%
#   filter(GEOID == "25025160101")
# 
# # isolate blockgroups within tract
# ma_blkgrp_1 <- ma_blkgrps %>%
#   filter(str_starts(GEOID, "25025160101"))
# 
# # tm_shape(ma_blkgrp_1) + tm_borders() +
# #   tm_shape(ma_blkgrp_1) + tm_text(text = "NAME")
# 
# # isolate a single block group within the tract
# ma_blkgrp_4 <- ma_blkgrp_1 %>%
#   filter(NAME == "1")
# 
# # create separate maps for each level
# # tmap_mode("plot")
# # first, create basemaps for each
# a_bm <- ma_counties %>%
#   st_transform(., crs = 4326) %>%
#   st_bbox(.) %>%
#   tmaptools::read_osm(., type = "esri")
# 
# # create a new bbox so that Chelsea map fills space
# bbox_new <- st_bbox(ma_cosub)
# 
# xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
# yrange <- bbox_new$ymax - bbox_new$ymin # range of y values
# 
# # bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
# # bbox_new[3] <- bbox_new[3] + (0.5 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.13 * yrange) # ymin - bottom
# bbox_new[4] <- bbox_new[4] + (0.1 * yrange) # ymax - top
# 
# bbox_new <- bbox_new %>%  # take the bounding box ...
#   st_as_sfc() # ... and make it a sf polygon
# 
# b_bm <- bbox_new %>%
#   st_transform(., crs = 4326) %>%
#   st_bbox(.) %>%
#   tmaptools::read_osm(., type = "esri")
# 
# # create a new bbox so that Block Group map fills space
# bbox_new2 <- st_bbox(ma_tracts_1)
# 
# xrange <- bbox_new2$xmax - bbox_new2$xmin # range of x values
# yrange <- bbox_new2$ymax - bbox_new2$ymin # range of y values
# 
# bbox_new2[1] <- bbox_new2[1] - (0.06 * xrange) # xmin - left
# bbox_new2[3] <- bbox_new2[3] + (0.06 * xrange) # xmax - right
# bbox_new2[2] <- bbox_new2[2] - (0.04 * yrange) # ymin - bottom
# bbox_new2[4] <- bbox_new2[4] + (0.01 * yrange) # ymax - top
# 
# bbox_new2 <- bbox_new2 %>%  # take the bounding box ...
#   st_as_sfc() # ... and make it a sf polygon
# 
# c_bm <- bbox_new2 %>%
#   st_transform(., crs = 4326) %>%
#   st_bbox(.) %>%
#   tmaptools::read_osm(., type = "esri")
# 
# # make the maps
# a <- tm_shape(a_bm) + tm_rgb() +
#   # tm_shape(ma_counties) + tm_fill(col = "grey88", alpha = 0.5) +
#   tm_shape(ma_cosub) + tm_fill(col = "#7570b3", alpha = 0.3) +
#   tm_shape(ma_cosub) + tm_borders(col = "#7570b3", lwd = 3) +
#   # tm_shape(ma_tracts_1) + tm_fill(col = "#d95f02") +
#   # tm_shape(ma_blkgrp_4) + tm_fill(col = "#1b9e77") +
#   tm_scale_bar(breaks = c(0,5,10), position = c("center","BOTTOM")) +
#   tm_layout(title = "Municipality", title.size = 0.9,
#             inner.margins = c(0.03,0.02,0.18,0.02))
# 
# b <- tm_shape(b_bm) + tm_rgb() +
#   tm_shape(ma_cosub, unit = "km") + tm_fill(col = "#7570b3",
#                                             alpha = 0.3) +
#   tm_shape(ma_tracts_1) + tm_fill(col = "#d95f02") +
#   # tm_shape(ma_blkgrp_4) + tm_fill(col = "#1b9e77") +
#   tm_shape(ma_cosub) + tm_borders(col = "#7570b3", lwd = 3) +
#   tm_credits("Tract", position = c(0.49,0.22)) +
#   tm_credits("Municipality", position = c(0.38,0.45),
#              col = "#7570b3", alpha = 0.7) +
#   tm_scale_bar(breaks = c(0,1,2), position = c("center","BOTTOM")) +
#   tm_layout(title = "Census Tract within\nMunicipality", title.size = 0.9,
#             inner.margins = c(0.01,0.02,0.16,0.02))
# 
# c <- tm_shape(c_bm, unit = "m", bbox = ma_tracts_1) + tm_rgb() +
#   tm_shape(ma_tracts_1, unit = "m") + tm_fill(col = "#d95f02",
#                                               alpha = 0.3) +
#   tm_shape(ma_blkgrp_4) + tm_fill(col = "#1b9e77") +
#   tm_shape(ma_tracts_1) + tm_borders(col = "#d95f02", lwd = 3) +
#   tm_credits("Block\nGroup", position = c(0.19,0.41)) +
#   tm_credits("Census Tract", position = c(0.38,0.30),
#              col = "#d95f02", alpha = 0.7) +
#   tm_scale_bar(breaks = c(0,250,500), position = c("center","BOTTOM")) +
#   tm_layout(title = "Census Block Group\nwithin Census Tract", title.size = 0.9,
#             inner.margins = c(0.1,0,0.21,0))
# 
# tmap_arrange(a,b,c, ncol = 3) %>%
#   tmap_save(., filename = "Images/pub/census_geog.jpeg",
#             width = 6, height = 3, units = "in", dpi = 1000)

knitr::include_graphics("Images/pub/census_geog.jpeg")
```


```{r MAUPtab, echo=FALSE, message=FALSE, warning=FALSE, comment="", fig.align='center'}
library(tidycensus)
library(tidyverse)
library(sf)
library(kableExtra)
# BGstats <- get_acs(geography = "block group", 
#                                variables = "B01001_001", 
#                                state = "MA", geometry = TRUE,
#                    show_call = FALSE) %>% 
#   mutate(sqkm = as.numeric(st_area(.))/10^6) %>% 
#   as.data.frame() %>% 
#   select(estimate, sqkm) %>% 
#   filter(sqkm > 0 & estimate > 0) %>% 
#   summarize(count = n(), minArea = min(sqkm), medArea = median(sqkm),
#             avgArea = mean(sqkm), maxArea = max(sqkm), 
#             minPop = min(estimate), medPop = median(estimate),
#             avgPop = mean(estimate), maxPop = max(estimate))
# 
# TRstats <- get_acs(geography = "tract", 
#                                variables = "B01001_001", 
#                                state = "MA", geometry = TRUE,
#                    show_call = FALSE) %>% 
#   mutate(sqkm = as.numeric(st_area(.))/10^6) %>% 
#   as.data.frame() %>% 
#   select(estimate, sqkm) %>% 
#   filter(sqkm > 0 & estimate > 0) %>% 
#   summarize(count = n(), minArea = min(sqkm), medArea = median(sqkm),
#             avgArea = mean(sqkm), maxArea = max(sqkm), 
#             minPop = min(estimate), medPop = median(estimate),
#             avgPop = mean(estimate), maxPop = max(estimate))
# 
# COstats <- get_acs(geography = "county subdivision", 
#                                variables = "B01001_001", 
#                                state = "MA", geometry = TRUE,
#                    show_call = FALSE) %>% 
#   mutate(sqkm = as.numeric(st_area(.))/10^6) %>% 
#   as.data.frame() %>% 
#   select(estimate, sqkm) %>% 
#   filter(sqkm > 0 & estimate > 0) %>% 
#   summarize(count = n(), minArea = min(sqkm), medArea = median(sqkm),
#             avgArea = mean(sqkm), maxArea = max(sqkm), 
#             minPop = min(estimate), medPop = median(estimate),
#             avgPop = mean(estimate), maxPop = max(estimate))
# 
# # compile into one df and format and save
# geoStats <- rbind(BGstats,TRstats,COstats) %>% 
#   cbind(data.frame(Geography = c("Block Group", "Tract",
#                                  "Municipality"))) %>% 
#   select(Geography, everything()) %>% 
#   mutate(minArea = round(minArea, 2), medArea = round(medArea,1), avgArea = round(avgArea,1), maxArea = round(maxArea,1), avgPop = round(avgPop,0), medPop = round(medPop,0))
# 
# save(geoStats, file = "Data/geoStats.Rds")

load("Data/geoStats.Rds")

geoStats %>% kable(longtable = T, booktabs = T, escape = F,
                    format.args = list(big.mark = ','),
                    caption = "Comparison of areas and populations for 
        units of analysis in Massachusetts", align = "r", 
        col.names = c("Geography", "count", "min", "med", "avg", 
                      "max", "min", "med", "avg", "max")) %>% 
  add_header_above(c(" " = 2, "Area (km$\\^{2}$)" = 4, "Population" = 4), escape = F) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```


# Results
```{r data, include=FALSE}
# load packages and data
library(tidyverse)
library(sf)
library(tigris)
library(tmap)
library(tmaptools)
library(ggplot2)
library(foreign) # for reading in dbf
# library(tidytext) # for reorder_within and scale_x_reorder in ggplot2 facets
library(kableExtra)
# library(sp)
# library(spdep)
# library(bibtex)

# Load demographic and gas leaks data
load("Data/Demographics.rds")
load("Data/HEET2019Leaksv2.rds")
ma_blkgrps <- readRDS("Data/ma_blkgrps2019.Rds")
ma_tracts <- readRDS("Data/ma_tracts2019.Rds")
ma_cosub <- readRDS("Data/ma_cosub2019.Rds")

ppLeakDensity_df_bg <- readRDS("Data/ppLeakDensity_df_blkgrps2019.Rds")
ppLeakDensity_df_tr <- readRDS("Data/ppLeakDensity_df_tracts2019.Rds")
ppLeakDensity_df_co <- readRDS("Data/ppLeakDensity_df_cosubs2019.Rds")

ppLeakDensityJoinedU_bg <- readRDS("Data/ppLeakDensityJoinedU_BG2019.Rds")
ppLeakDensityJoinedU_tr <- readRDS("Data/ppLeakDensityJoinedU_TR2019.Rds")
ppLeakDensityJoinedU_co <- readRDS("Data/ppLeakDensityJoinedU_CO2019.Rds")

# set common CRS
unrepaired2019 <- st_transform(unrepaired2019, crs = st_crs(ma_blkgrps))
repaired2019 <- st_transform(repaired2019, crs = st_crs(ma_blkgrps))

# Set up recurring map elements
# Load natural gas utility service areas dbf from MassGIS and join to MassGIS towns layer
ng_dbf <- read.dbf("Data/pubutil/TOWNS_POLY_UTILITIES.dbf")

ng_service_areas <- st_read(dsn = "Data/townssurvey_shp",
                            layer = "TOWNSSURVEY_POLYM") %>% 
  select(-TOWN) %>% 
  left_join(., ng_dbf, by = "TOWN_ID") %>% 
  select(TOWN, GAS, GAS_LABEL) %>%  
  st_transform(., crs = st_crs(ma_blkgrps)) %>% 
  st_make_valid()

# isolate No gas areas and municipal
ng_nogas_muni <- ng_service_areas %>% 
  filter(GAS_LABEL %in% c("No Natural Gas Service","Municipal")) %>% 
  group_by(GAS_LABEL) %>%
  summarize() %>% 
  mutate(GAS_LABEL = as.character(GAS_LABEL))

# consolidate into basic gas service districts
ng_service_areas2 <- ng_service_areas %>% 
  filter(!GAS_LABEL %in% c("No Natural Gas Service","Municipal")) %>% 
  group_by(GAS_LABEL) %>% 
  summarize(TownCnt = n()) %>% 
  mutate(ID = case_when(
    GAS_LABEL == "National Grid" ~ "NG",
    GAS_LABEL == "Blackstone Gas Company" ~ "BGC",
    GAS_LABEL == "Columbia Gas" ~ "CG",
    GAS_LABEL == "Eversource Energy" ~ "EV",
    GAS_LABEL == "Columbia Gas, Eversource Energy" ~ "CG,EV",
    GAS_LABEL == "Unitil" ~ "UN",
    GAS_LABEL == "National Grid, Unitil" ~ "NG,UN",
    GAS_LABEL == "Eversource Energy, National Grid" ~ "EV,NG",
    GAS_LABEL == "The Berkshire Gas Company" ~ "BG",
    GAS_LABEL == "Columbia Gas, Blackstone Gas Company" ~ "CG,BGC",
    GAS_LABEL == "Liberty Utilities" ~ "LU",
    GAS_LABEL == "Colonial Gas" ~ "NG",
    GAS_LABEL == "Columbia Gas, National Grid" ~ "CG,NG"
  ))

# join ID to original municipalities within ng_service_area to use centroids as labels; restrict to central municipalities within territories
ng_service_labels <- ng_service_areas2 %>% 
  as.data.frame() %>% 
  select(ID, GAS_LABEL) %>% 
  inner_join(ng_service_areas, ., by = "GAS_LABEL") %>% 
  filter(TOWN %in% c("LANESBOROUGH","DEERFIELD","NORTHAMPTON","LUDLOW",
                     "NORTH BROOKFIELD","OXFORD",
                     "WESTMINSTER","CARLISLE","ANDOVER","WESTBOROUGH",
                     "BLACKSTONE","MILTON",
                     "NORWOOD","WEST BRIDGEWATER","PLAINVILLE","PLYMOUTH",
                     "WESTPORT","ACUSHNET","WAREHAM","BARNSTABLE",
                     "SWANSEA","MARSHFIELD","COHASSET",
                     "HAMILTON","WEST NEWBURY","DOVER","DUNSTABLE")) %>% 
  st_centroid(., of_largest_polygon = TRUE)

ng_service_labels2 <- ng_service_areas2 %>% 
  as.data.frame() %>% 
  select(ID, GAS_LABEL) %>% 
  inner_join(ng_service_areas, ., by = "GAS_LABEL") %>% 
  filter(TOWN %in% c("LEICESTER","LUNENBURG","MENDON","BELLINGHAM","HANSON",
                     "WAYLAND","NATICK","BOSTON","SOMERVILLE")) %>% 
  st_centroid(., of_largest_polygon = TRUE)

ng_service_labels3 <- ng_service_areas2 %>% 
  as.data.frame() %>% 
  select(ID, GAS_LABEL) %>% 
  inner_join(ng_service_areas, ., by = "GAS_LABEL") %>% 
  filter(TOWN == "BOSTON") %>% 
  st_centroid(., of_largest_polygon = TRUE)

# separate MA for cropping
ma_state_sf <- states(cb = TRUE) %>%
  filter(STUSPS == "MA")

# grab municipal boundaries
# ma_towns_sf <- county_subdivisions(state = "MA", cb = TRUE) %>% 
#   st_transform(., crs = 26986) # breaking periodically

# Download file directly from https://www2.census.gov/geo/tiger/TIGER2019/
ma_towns_sf <- st_read(dsn = "Data/tl_2019_25_cousub", 
                       layer = "tl_2019_25_cousub") %>% 
  st_transform(., crs = 26986)

# create point layer of towns for context
ma_towns_sf_pts <- ma_towns_sf %>% 
  # county_subdivisions(state = "MA", cb = TRUE) %>% 
  filter(NAME %in% c("Boston",
                     "Lawrence",
                     "Lowell",
                     "Brockton",
                     "Worcester",
                     "Springfield",
                     "Pittsfield", 
                     "Stockbridge",
                     "Fall River",
                     "West Yarmouth",
                     "Lynn",
                     "Randolph",
                     "Webster",
                     "Attleboro",
                     "Medford",
                     "Amherst",
                     "Quincy",
                     "Weymouth Town",
                     "Nantucket")) %>% 
  mutate(NAME = recode(NAME, "Weymouth Town" = "Weymouth")) %>% 
  st_transform(., crs = 26986) %>% 
  st_centroid(of_largest_polygon = TRUE)

# create a separate point for Newton so that it can be repositioned
newton <- ma_towns_sf %>% 
  # county_subdivisions(state = "MA", cb = TRUE) %>%  
  filter(NAME %in% c("Fitchburg","Newton","Edgartown","New Bedford")) %>% 
  st_transform(., crs = 26986) %>% 
  st_centroid(of_largest_polygon = TRUE)
```
In calendar year 2019, there were `r format(nrow(unrepaired2019)+nrow(repaired2019), big.mark = ",")` natural gas leaks in the distribution system reported by six of the seven investor-owned utilities in Massachusetts that filed leak reports with the DPU. These utilities represent 98% of natural gas customers in the state [@AGA] (see Figure \@ref(fig:mapUtilities)).[^d] 

```{r mapUtilities, fig.align='center', fig.cap= "Massachusetts Natural Gas Utility Territories, 2019", cache=TRUE}
# Map of utility territories

# create two subsets of territories to break up legend
# first extract names of utilities and convert to sorted text
util_names1 <- unique(ng_service_areas$GAS_LABEL) %>% 
  as.character(.) %>% 
  recode(., "Colonial Gas" = "National Grid") %>% 
  unique(.) %>% 
  sort(.) %>% 
  .[1:7]

util_names2 <- unique(ng_service_areas$GAS_LABEL) %>% 
  as.character(.) %>% 
  recode(., "Colonial Gas" = "National Grid") %>%
  unique(.) %>% 
  sort(.) %>% 
  .[8:14]
# create subsets of sf utility territories
ng_service_areas_BE <- ng_service_areas %>% 
  filter(GAS_LABEL %in% util_names1) %>% 
  mutate(GAS_LABEL = as.character(GAS_LABEL))
# fold Colonial Gas into National Grid
ng_service_areas_LU <- ng_service_areas %>% 
  mutate(GAS_LABEL = if_else(as.character(GAS_LABEL) == "Colonial Gas", "National Grid", as.character(GAS_LABEL))) %>% 
  filter(GAS_LABEL %in% util_names2)%>% 
  mutate(GAS_LABEL = as.character(GAS_LABEL))

# set palette of colors for all territories based on full sf
# "#1f78b4" was Colonial Gas, but folded into National Grid
myPal <- c("#a6cee3", "#b2df8a", "#33a02c", "#fb9a99",
           "#e31a1c", "#fdbf6f", "#ff7f00", "blue", "ivory2", "#cab2d6",
           "#6a3d9a", "gray88", "#ffff99", "#b15928")

# create palette subsets for legend only layouts
myPal1 <- myPal[1:7]
myPal2 <- myPal[8:14]

# create main map object of utility territories
m_Util <- tm_shape(ng_service_areas_BE, bbox = ng_service_areas2) + 
  tm_fill(col = "GAS_LABEL", palette = myPal1, title = "") +
  tm_shape(ng_service_areas_LU) + 
  tm_fill(col = "GAS_LABEL", palette = myPal2, title = "") + 
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("center","BOTTOM")) +
  tm_layout(legend.position = c(0.02,0.04), 
            legend.stack = "horizontal", legend.text.size = 0.5,
            legend.width = 0.65, title = "Massachusetts Natural Gas Utility Territories, 2019",
            title.size = 0.9,
            inner.margins = c(0.02,0.02,0.09,0.02))


# create inset map of US states for context
us_states <- tigris::states(cb = TRUE) %>% 
  filter(!NAME %in% c("Alaska", "Hawaii", "Puerto Rico", "American Samoa", "Commonwealth of the Northern Mariana Islands", "United States Virgin Islands"))
ma_state <- us_states %>% 
  filter(NAME == "Massachusetts")

# create modified bbox to focus inset on contiguous US
bbox_new <- st_bbox(us_states)
bbox_new["xmax"] <- -67
bbox_new["ymin"] <- 25

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

# create the inset map object
inset <- tm_shape(us_states, bbox = bbox_new) + 
  tm_polygons(col = "grey88", border.col = "grey", alpha = 0.7) +
  tm_shape(ma_state) + tm_fill(col = "red")

# create viewport to specify dimensions of inset
vp_inset <- grid::viewport(x = 0.795, y = 0.85, width = 0.5, 
                           height = 0.25)

tmap_save(m_Util, filename = "Images/pub/m_Util.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = inset,
          insets_vp = vp_inset)

knitr::include_graphics("Images/pub/m_Util.jpeg")
```

Repairs were reported for `r format(nrow(repaired2019), big.mark = ",")` (`r round(nrow(repaired2019)/(nrow(unrepaired2019)+nrow(repaired2019)),2)*100`%) of these, leaving `r format(nrow(unrepaired2019), big.mark = ",")` (`r round(nrow(unrepaired2019)/(nrow(unrepaired2019)+nrow(repaired2019)),2)*100`%) unrepaired leaks across the state as of December 31, 2019 (see Table  \@ref(tab:tabLeaks)). 

[^d]: No gas leak data was available from the DPU for Blackstone Gas Company.

```{r tabLeaks, echo=FALSE, message=FALSE, warning=FALSE, fig.align='center'}
# table of leak stats by class and repair status
unrepaired_Sumdf <- unrepaired2019 %>% as.data.frame() %>% 
  group_by(Class) %>% summarize(Unrepaired = n())

repaired_Sumdf <- repaired2019 %>% as.data.frame() %>% 
  group_by(Class) %>% summarize(Repaired = n())

total_Sumdf <- left_join(unrepaired_Sumdf, repaired_Sumdf, by = "Class") %>% 
  drop_na() %>%
  mutate(Total = Unrepaired + Repaired) %>% 
  bind_rows(summarize(.,
                      across(where(is.numeric), sum),
                      across(where(is.character), ~ "Total"))) %>% 
  mutate(PctUnrepaired = Unrepaired/Total*100,
         PctRepaired = Repaired/Total*100,
         PctTotal = Total/Total[Class == "Total"]*100) %>%
  select(Class, Unrepaired, PctUnrepaired, Repaired, PctRepaired, 
         Total, PctTotal)

total_Sumdf %>% 
  mutate(PctUnrepaired = paste0(round(PctUnrepaired,1),"%"),
         PctRepaired = paste0(round(PctRepaired,1),"%"),
         PctTotal = paste0(round(PctTotal,1),"%")) %>% 
  kable(., longtable = T, booktabs = T,
        format.args = list(big.mark = ','), 
        caption = "Gas leaks by class and repair status, 2019", align = "r", 
        # digits = c(0,0,1,0,1,0,1), 
        col.names = c("Class","Count","Percent","Count","Percent",
                      "Total","Pct Total")) %>% 
  add_header_above(., c(" ", "Unrepaired Leaks" = 2, "Repaired Leaks" = 2, 
                        " " = 2)) %>% 
  row_spec(4, bold = T) %>% 
  kable_styling(latex_options = c("repeat_header"))
```

Reported gas leaks were nearly ubiquitous across the state (see Figures \@ref(fig:mapAllLeaks) - \@ref(fig:mapAllLeaksCO)). At least one gas leak was reported in `r round((ma_cosub %>% filter(AllLeaks2019>0) %>% nrow())/nrow(ma_cosub)*100,0)`% of the municipalities within natural gas service territories, and `r round((ma_cosub %>% filter(unrepaired2019total>0) %>% nrow())/nrow(ma_cosub)*100,0)`% had at least one unrepaired gas leak at the end of 2019. Over `r round((ma_blkgrps %>% filter(AllLeaks2019>0) %>% nrow())/nrow(ma_blkgrps)*100,0)`% of Census Block Groups had at least one known gas leak, and approximately `r round((ma_blkgrps %>% filter(unrepaired2019total>0) %>% nrow())/nrow(ma_blkgrps)*100,0)`% had at least one unrepaired gas leak at the end of 2019. However, there was significant variation in geographic distribution and concentration of these leaks.

```{r mapAllLeaks, fig.align='center', fig.cap= "Massachusetts Gas Leaks and Utility Territories, 2019, aggregated by 1km hexagon tessellations", cache=TRUE}
# Map of leak counts by hexagons and utility territories

# create a hexagonal grid and create index column
gridCnt <- st_make_grid(x = ma_blkgrps, cellsize = 1000, square = FALSE) %>% 
  st_as_sf() %>% 
  mutate(index = row_number())

# spatially join to grid ids to leaks, sum aggregate numbers of leaks per grid index, and then join sums back to hexagons for mapping and analysis
gridCnt <- unrepaired2019 %>% 
  st_join(., gridCnt) %>% 
  st_drop_geometry() %>% 
  group_by(index) %>% 
  summarize(unrepaired = n()) %>% 
  left_join(gridCnt, ., by = "index") %>% 
  replace_na(list(unrepaired = 0))

gridCnt <- repaired2019 %>% 
  st_join(., gridCnt) %>% 
  st_drop_geometry() %>% 
  group_by(index) %>% 
  summarize(repaired = n()) %>% 
  left_join(gridCnt, ., by = "index") %>% 
  replace_na(list(repaired = 0))

# create column with total leak points and clip to MA
gridCnt <- gridCnt %>% 
  mutate(total = unrepaired + repaired) %>% 
  crop_shape(., ma_blkgrps, polygon = TRUE) %>% 
  st_make_valid()

gridCntT <- gridCnt %>% 
  filter(total > 0)

m_gridTotal <- tm_shape(gridCntT, unit = "km", bbox = ma_state_sf) + 
  tm_fill(col = "total", palette = "YlOrRd",
          style = "fisher",
          # breaks = c(1,5,10,20,40,80), 
          legend.format = list(digits = 0),
          title = "Number\nof leaks") +
  tm_shape(ng_nogas_muni, bbox = ma_state_sf) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.055,0.1),
            title = "Massachusetts Gas Leaks and Utility Territories, 2019",
            title.size = 0.9,
            inner.margins = c(0.02,0.02,0.09,0.02)
  ) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m_gridTotal, filename = "Images/pub/m_gridTotalUtil.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_gridTotalUtil.jpeg")
```

```{r mapAllLeaksBG, fig.align='center', fig.cap= "Massachusetts gas leaks per square kilometer by Census Block Group, 2019", cache=TRUE}
# create map of leak densities by BG
blkgrpsLeaksAll <- ma_blkgrps %>% 
  filter(AllLeaks2019_sqkm > 0)

m1 <- tm_shape(blkgrpsLeaksAll, bbox = ma_state_sf) + 
  tm_fill(col = "AllLeaks2019_sqkm", style = "fisher", palette = "YlOrRd", 
          title = "Leaks per SqKm", legend.format = list(digits = 2)) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_blkgrps) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + 
  tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.15),
            title = "Massachusetts Gas Leak Density by Census Block Group and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m1, filename = "Images/pub/m_blkgrpTotal.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_blkgrpTotal.jpeg")
```

```{r mapAllLeaksTR, fig.align='center', fig.cap= "Massachusetts gas leaks per square kilometer by Census Tract, 2019",cache=TRUE}
tractsLeaksAll <- ma_tracts %>% 
  filter(AllLeaks2019_sqkm > 0)

m2 <- tm_shape(tractsLeaksAll, bbox = ma_state_sf) + 
  tm_fill(col = "AllLeaks2019_sqkm", 
          style = "fisher",
          palette = "YlOrRd", 
          title = "Leaks per SqKm",
          legend.format = list(digits = 2)) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_tracts) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.15),
            title = "Massachusetts Gas Leak Density by Census Tract and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m2, filename = "Images/pub/m_tractTotal.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_tractTotal.jpeg")
```

```{r mapAllLeaksCO, fig.align='center', fig.cap= "Massachusetts gas leaks per square kilometer by municipality, 2019",cache=TRUE}
cosubLeaksAll <- ma_cosub %>% 
  filter(AllLeaks2019_sqkm > 0)

m3 <- tm_shape(cosubLeaksAll, bbox = ma_state_sf) + 
  tm_fill(col = "AllLeaks2019_sqkm", style = "fisher", palette = "YlOrRd", 
          title = "Leaks per SqKm",legend.format = list(digits = 2)) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_cosub) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.15),
            title = "Massachusetts Gas Leak Density by Municipality and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m3, filename = "Images/pub/m_cosubTotal.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_cosubTotal.jpeg")
```

```{r tabBGcount, fig.align='center', include=FALSE}
# table of leak count stats by block group
# list of stats to compute across
summary_stats <- list(Min = ~min(., na.rm = T),
                      Med = ~median(., na.rm = T),
                      Avg = ~mean(., na.rm = T),
                      Max = ~max(., na.rm = T))
# table of stats by block group
unrep_blkgrp_cnt <- ma_blkgrps %>% 
  as.data.frame() %>% 
  select(starts_with("unrepaired")) %>% 
  pivot_longer(cols = everything(), names_to = "Class", values_to = "Count") %>% 
  group_by(Class) %>% 
  summarize(across(.cols = Count, summary_stats, 
                             .names = "{.fn}U")) %>% 
  mutate(Class = recode(Class, "unrepaired2019total" = "All",
                        "unrepaired2019totalC1" = "1",
                        "unrepaired2019totalC2" = "2",
                        "unrepaired2019totalC3" = "3")) %>% 
  arrange(Class)

rep_blkgrp_cnt <- ma_blkgrps %>% 
  as.data.frame() %>% 
  select(starts_with("repaired")) %>% 
  pivot_longer(cols = everything(), names_to = "Class", values_to = "Count") %>% 
  group_by(Class) %>% 
  summarize(across(.cols = Count, summary_stats, 
                   .names = "{.fn}R")) %>% 
  mutate(Class = recode(Class, "repaired2019total" = "All",
                        "repaired2019totalC1" = "1",
                        "repaired2019totalC2" = "2",
                        "repaired2019totalC3" = "3")) %>% 
  arrange(Class)

all_blkgrp_cnt <- ma_blkgrps %>% 
  as.data.frame() %>% 
  select(AllLeaks2019:AllLeaks2019C3) %>% 
  pivot_longer(cols = everything(), names_to = "Class", values_to = "Count") %>% 
  group_by(Class) %>% 
  summarize(across(.cols = Count, summary_stats, 
                   .names = "{.fn}T")) %>% 
  mutate(Class = recode(Class, "AllLeaks2019" = "All",
                        "AllLeaks2019C1" = "1",
                        "AllLeaks2019C2" = "2",
                        "AllLeaks2019C3" = "3")) %>% 
  arrange(Class)

# join together
list(unrep_blkgrp_cnt,rep_blkgrp_cnt,all_blkgrp_cnt) %>% 
  reduce(., left_join, by = "Class") %>% 
  kable(., longtable = T, booktabs = T,
        # format.args = list(big.mark = ','), 
        caption = "Gas leak counts by class per Census Block Group, 2019", 
        align = "r", 
        digits = c(0,0,0,1,0,0,0,1,0,0,0,1),
        col.names = c("Class","Min","Med","Avg","Max",
                      "Min","Med","Avg","Max",
                      "Min","Med","Avg","Max")) %>% 
  add_header_above(., c(" ", "Unrepaired Leaks" = 4, "Repaired Leaks" = 4, 
                        "Total Leaks" = 4)) %>% 
  # row_spec(4, bold = T) %>% 
  kable_styling(latex_options = c("repeat_header"))
```

```{r calcUtility, include=FALSE}
# create df of counts per utility
unrepairedUSum_df <- unrepaired2019 %>% 
  as.data.frame() %>% 
  mutate(Utility = recode(Utility, 
                          "Fitchburg" = "Unitil/Fitchburg Gas",
                          "Eversource" = "Eversource Energy")) %>% 
  group_by(Utility) %>% 
  summarize(Unrepaired = n())

repairedUSum_df <- repaired2019 %>% 
  as.data.frame() %>% 
  mutate(Utility = recode(Utility, 
                          "Fitchburg" = "Unitil/Fitchburg Gas",
                          "Eversource" = "Eversource Energy")) %>% 
  group_by(Utility) %>% 
  summarize(Repaired = n())

totalUSum_df <- left_join(unrepairedUSum_df, repairedUSum_df, by = "Utility") %>% 
  mutate(Total = Unrepaired + Repaired) %>% 
  bind_rows(summarize(.,
                      across(where(is.numeric), sum),
                      across(where(is.character), ~ "Total"))) %>% 
  mutate(PctUnrepaired = Unrepaired/Total*100,
         PctRepaired = Repaired/Total*100,
         PctTotal = Total/Total[Utility == "Total"]*100) %>% 
  select(Utility, Unrepaired, PctUnrepaired, Repaired, PctRepaired, Total, PctTotal)
```

The vast majority of all leaks were located in the eastern half of the state, particularly in the greater Boston region, where most of the population resides, but also within major cities throughout the state. Moreover, most of these leaks (`r totalUSum_df %>% filter(Utility == "National Grid") %>% select(PctTotal) %>% pull() %>% round(.,0)`%) occurred within the service territory of National Grid, which is the largest natural gas utility in the state, serving Boston and surrounding communities (see Figure \@ref(fig:figLeaksUtility)).

```{r tabBGdensity, fig.align='center', include=FALSE}
# table of leak density stats by block group
# list of stats to compute across
summary_stats <- list(Min = ~min(., na.rm = T),
                      Med = ~median(., na.rm = T),
                      Avg = ~mean(., na.rm = T),
                      Max = ~max(., na.rm = T))
# table of stats by block group
unrep_blkgrp_dns <- ma_blkgrps %>% 
  as.data.frame() %>% 
  select(starts_with("leaks_sqkm")) %>% 
  pivot_longer(cols = everything(), names_to = "Class", values_to = "Count") %>% 
  group_by(Class) %>% 
  summarize(across(.cols = Count, summary_stats, 
                   .names = "{.fn}U")) %>% 
  mutate(Class = recode(Class, "leaks_sqkm" = "All",
                        "leaks_sqkmC1" = "1",
                        "leaks_sqkmC2" = "2",
                        "leaks_sqkmC3" = "3")) %>% 
  arrange(Class)

rep_blkgrp_dns <- ma_blkgrps %>% 
  as.data.frame() %>% 
  select(starts_with("REPleaks_sqkm")) %>% 
  pivot_longer(cols = everything(), names_to = "Class", values_to = "Count") %>% 
  group_by(Class) %>% 
  summarize(across(.cols = Count, summary_stats, 
                   .names = "{.fn}R")) %>% 
  mutate(Class = recode(Class, "REPleaks_sqkm" = "All",
                        "REPleaks_sqkmC1" = "1",
                        "REPleaks_sqkmC2" = "2",
                        "REPleaks_sqkmC3" = "3")) %>% 
  arrange(Class)

all_blkgrp_dns <- ma_blkgrps %>% 
  as.data.frame() %>% 
  select(AllLeaks2019_sqkm:AllLeaks2019C3_sqkm) %>% 
  pivot_longer(cols = everything(), names_to = "Class", values_to = "Count") %>% 
  group_by(Class) %>% 
  summarize(across(.cols = Count, summary_stats, 
                   .names = "{.fn}T")) %>% 
  mutate(Class = recode(Class, "AllLeaks2019_sqkm" = "All",
                        "AllLeaks2019C1_sqkm" = "1",
                        "AllLeaks2019C2_sqkm" = "2",
                        "AllLeaks2019C3_sqkm" = "3")) %>% 
  arrange(Class)

# join together
list(unrep_blkgrp_dns,rep_blkgrp_dns,all_blkgrp_dns) %>% 
  reduce(., left_join, by = "Class") %>% 
  kable(., longtable = T, booktabs = T,
        # format.args = list(big.mark = ','), 
        caption = "Gas leak density (per sqkm) by class per Census Block Group, 2019", 
        align = "r", 
        digits = 2,
        col.names = c("Class","Min","Med","Avg","Max",
                      "Min","Med","Avg","Max",
                      "Min","Med","Avg","Max")) %>% 
  add_header_above(., c(" ", "Unrepaired Leaks" = 4, "Repaired Leaks" = 4, 
                        "Total Leaks" = 4)) %>% 
  # row_spec(4, bold = T) %>% 
  kable_styling(latex_options = c("repeat_header"))
```

```{r figLeaksUtility, fig.align='center', fig.cap="Frequency of reported repaired and unrepaired leaks by natural gas utility in Massachusetts, 2019"}
# create a stacked bar chart of leaks by utility. shows significant difference in magnitude. 
unrepairedUSum_df <- unrepaired2019 %>% 
  as.data.frame() %>% 
  mutate(Utility = recode(Utility, 
                          "Fitchburg" = "Unitil/Fitchburg Gas",
                          "Eversource" = "Eversource Energy")) %>% 
  group_by(Utility) %>% 
  summarize(Unrepaired = n())

repairedUSum_df <- repaired2019 %>% 
  as.data.frame() %>% 
  mutate(Utility = recode(Utility, 
                          "Fitchburg" = "Unitil/Fitchburg Gas",
                          "Eversource" = "Eversource Energy")) %>% 
  group_by(Utility) %>% 
  summarize(Repaired = n())

totalUSum_df <- left_join(unrepairedUSum_df, repairedUSum_df, by = "Utility") %>% 
  mutate(Total = Unrepaired + Repaired) %>% 
  bind_rows(summarize(.,
                      across(where(is.numeric), sum),
                      across(where(is.character), ~ "Total"))) %>% 
  mutate(PctUnrepaired = Unrepaired/Total*100,
         PctRepaired = Repaired/Total*100,
         PctTotal = Total/Total[Utility == "Total"]*100) %>% 
  select(Utility, Unrepaired, PctUnrepaired, Repaired, PctRepaired, Total, PctTotal)

totalUSum_df %>% 
  filter(Utility != "Total") %>% 
  select(-c(Total, starts_with("Pct"))) %>% 
  pivot_longer(., cols = c(Unrepaired,Repaired), 
               names_to = "status", values_to = "leaks") %>% 
    ggplot(aes(x = reorder(Utility, leaks), y = leaks, 
             fill = status)) +
    geom_bar(position = "stack", stat = "identity") +
  scale_fill_manual("legend", values = c("Repaired" = "#7CB5EC", 
                                         "Unrepaired" = "#F7A35C")) + 
  coord_flip() +
  theme_minimal() +
  theme(legend.title = element_blank(), legend.position = "top") +
  scale_y_continuous(labels=function(x) format(x, big.mark = ",", scientific = FALSE)) +
  guides(fill = guide_legend(reverse = TRUE)) +
  labs(x = NULL, 
       y = "Reported Leaks",
       title = "Reported leaks by utility in 2019 across Massachusetts")

ggsave("Images/pub/LeaksbyUtility.jpeg", dpi = 1000)

# knitr::include_graphics("Images/pub/LeaksbyUtility.jpeg")
```


## Exposure
### Statewide exposure

Exposure to unrepaired leaks demonstrates consistently strong differences both geographically and demographically. White residents and people over 64 consistently experience below average population-weighted leak densities for all classes or grades of unrepaired leaks. This is the most invariant finding and is unaffected by the scale or unit of analysis, and true in every utility territory except Berkshire Gas in western Massachusetts. Conversely, People of Color (particularly Asians and Blacks), as well as limited-English speaking households, consistently experience higher densities of unrepaired leaks. Moreover, race and English-language ability exhibit more consistent leak exposure disparities than income indicators in every utility territory save for Berkshire Gas. 

Figure \@ref(fig:figDotUnrepaired) shows the relative exposure by leak class for unrepaired leaks for each population group at the Census Block Group (square), Census Tract (triangle), and municipality (circle) levels. For all leak classes, White people experience population-weighted mean relative exposures to unrepaired leaks that range from `r round(1-ppLeakDensity_df_co %>% filter(Group == "White") %>% select(wLeaksRR) %>% pull(),2)*100`% (RE `r round(ppLeakDensity_df_co %>% filter(Group == "White") %>% select(wLeaksRR) %>% pull(),2)` at municipality level) to `r round(1-ppLeakDensity_df_bg %>% filter(Group == "White") %>% select(wLeaksRR) %>% pull(),2)*100`% (RE `r round(ppLeakDensity_df_bg %>% filter(Group == "White") %>% select(wLeaksRR) %>% pull(),2)` at Block Group level) below that of the general population. By contrast, Black people experience relative exposures of `r abs(round(1-ppLeakDensity_df_co %>% filter(Group == "Black") %>% select(wLeaksRR) %>% pull(),2))*100`% (RE `r round(ppLeakDensity_df_co %>% filter(Group == "Black") %>% select(wLeaksRR) %>% pull(),2)` at municipality level) to `r abs(round(1-ppLeakDensity_df_bg %>% filter(Group == "Black") %>% select(wLeaksRR) %>% pull(),2))*100`% (RE `r round(ppLeakDensity_df_bg %>% filter(Group == "Black") %>% select(wLeaksRR) %>% pull(),2)` at Block Group level) above that of the general population. The differences are starker for potentially more hazardous Class 1 and Class 2 leaks. 
For Class 1 unrepaired leak densities, Whites have a relative exposure of approximately `r round(1-ppLeakDensity_df_co %>% filter(Group == "White") %>% select(wLeaksRRC1) %>% pull(),2)*100`% (RE `r round(ppLeakDensity_df_co %>% filter(Group == "White") %>% select(wLeaksRRC1) %>% pull(),2)`) below that of the general population at all scales of analysis. By comparison, Blacks have exposures of `r abs(round(1-ppLeakDensity_df_bg %>% filter(Group == "Black") %>% select(wLeaksRRC1) %>% pull(),2))*100`% (RE `r round(ppLeakDensity_df_bg %>% filter(Group == "Black") %>% select(wLeaksRRC1) %>% pull(),2)` at Block Group level) to `r abs(round(1-ppLeakDensity_df_co %>% filter(Group == "Black") %>% select(wLeaksRRC1) %>% pull(),2))*100`% (`r round(ppLeakDensity_df_co %>% filter(Group == "Black") %>% select(wLeaksRRC1) %>% pull(),2)` at municipality level) above that of the general population. For Class 2 unrepaired leak densities, Whites have a relative exposure of `r round(1-ppLeakDensity_df_co %>% filter(Group == "White") %>% select(wLeaksRRC2) %>% pull(),2)*100`% (RE `r round(ppLeakDensity_df_co %>% filter(Group == "White") %>% select(wLeaksRRC2) %>% pull(),2)` at municipality level) to `r round(1-ppLeakDensity_df_tr %>% filter(Group == "White") %>% select(wLeaksRRC2) %>% pull(),2)*100`% (`r round(ppLeakDensity_df_tr %>% filter(Group == "White") %>% select(wLeaksRRC2) %>% pull(),2)` at Tract level) below that of the general population. Blacks have exposures of `r abs(round(1-ppLeakDensity_df_co %>% filter(Group == "Black") %>% select(wLeaksRRC2) %>% pull(),2))*100`% (RE `r round(ppLeakDensity_df_co %>% filter(Group == "Black") %>% select(wLeaksRRC2) %>% pull(),2)` at municipality level) to `r abs(round(1-ppLeakDensity_df_tr %>% filter(Group == "Black") %>% select(wLeaksRRC2) %>% pull(),2))*100`% (`r round(ppLeakDensity_df_tr %>% filter(Group == "Black") %>% select(wLeaksRRC2) %>% pull(),2)` at Tract level) above that of the general population. Residents in state-designated, limited English environmental justice communities (MA Limited English HH) exhibit the greatest relative exposures for higher hazard leaks; over `r abs(round(1-ppLeakDensity_df_bg %>% filter(Group == "MA_ENGLISH") %>% select(wLeaksRRC2) %>% pull(),0))*100`% above that of the general population.  

Differences in relative exposure vary only slightly by the unit or scale of analysis. In general, the differences in relative exposure are most pronounced at the Census Block Group level (smallest geographic unit of analysis) and are lowest at the municipality level (largest geographic unit of analysis). For Class 2 leaks, however, Block Group-level differences fall in between those of the Tract and municipality level analyses. Relative exposures for Class 1 and Class 2 leaks exhibit the largest spread across the scales of analysis. However, Class 3 leaks, which account for the majority of unrepaired leaks, show considerable overlap or agreement between the scales of analysis. 

```{r figDotUnrepaired, fig.align='center', fig.cap="Relative exposures to population-weighted mean unrepaired gas leak density at Census Block Group, Census Tract, and municipality scales"}
# dot graph of relative exposure for unrepaired leak density around zero line
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

# create a consistent factor order to keep groups in a consistent order across all graphs
group_orderBG <- ppLeakDensity_df_bg %>% 
  arrange(desc(wLeaksRR)) %>% 
  select(Group) %>% 
  pull()

group_orderTR <- ppLeakDensity_df_tr %>% 
  arrange(wLeaksRR) %>% 
  select(Group) %>% 
  pull()

group_orderCO <- ppLeakDensity_df_co %>% 
  arrange(wLeaksRR) %>% 
  select(Group) %>% 
  pull()

# Set up for dot graph to show all three scales
ppLeakDensity_df_bg2 <- ppLeakDensity_df_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>%
  pivot_longer(wLeaksRR:wLeaksRRC3, 
               names_to = "leakClass", values_to = "leakDensityBG") %>% 
  # filter(!Group %in% c("MA_ENGLISH","MA_MINORITY21","MA_INCOME21")) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income")) %>%
  mutate(leakClass = recode(leakClass, "wLeaksRR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeaksRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeaksRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeaksRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, leakClass),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, leakClass, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensity_df_tr2 <- ppLeakDensity_df_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>%
  pivot_longer(wLeaksRR:wLeaksRRC3, 
               names_to = "leakClass", values_to = "leakDensityTR") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wLeaksRR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeaksRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeaksRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeaksRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, leakClass),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, leakClass, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensity_df_co2 <- ppLeakDensity_df_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(wLeaksRR:wLeaksRRC3, 
               names_to = "leakClass", values_to = "leakDensityCO") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wLeaksRR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeaksRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeaksRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeaksRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, leakClass),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, leakClass, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensity_df_bg2 <- ppLeakDensity_df_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensity_df_bg2,.)

# join them all together
ppLeakDensity_df_all <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                             ppLeakDensity_df_co2) %>% 
  lapply(., function(x) mutate(x, ClassGroup = paste(leakClass,Group))) %>%
  reduce(., full_join, by = "ClassGroup") %>% 
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensity_df_all2 <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                              ppLeakDensity_df_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
ppLeakDensity_df_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensity_df_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ leakClass, scales = "free_y") +
  labs(x = NULL, 
       y = expression(paste("Ratio of group population-weighted mean leak density (leaks/", 
                            km^2, ")", " to total population-weighted mean",sep = "")),
       title = "Relative Exposure to Unrepaired Gas Leaks in 2019 across Massachusetts", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotGraphallscalesRE.jpeg", dpi = 1000)
```


### Utility-specific exposure

When broken out by utility, relative exposure to population-weighted unrepaired leak density is largely consistent with the statewide pattern, with a few notable variations (see Figure \@ref(fig:figDotUnrepairedUtility)). For five of the six utilities examined here, Whites consistently show the lowest relative exposure, while limited English speaking households, People of Color (especially Blacks), and lower income groups (including renters) are more exposed than the general population. These differences are greatest for Columbia Gas and Unitil/Fitchburg Gas. The specific groups with greatest relative exposures vary between utility territories. For National Grid and Eversource Energy, Blacks show the highest relative exposures, for Columbia Gas it is limited English speaking households and Hispanics, and for Berkshire Gas and Unitil it is residents living in MA Low Income environmental justice communities.

```{r figDotUnrepairedUtility, fig.align='center', fig.cap="Relative exposures to population-weighted mean unrepaired gas leak density by utility at Census Block Group, Census Tract, and municipality scales"}
# Set up for dot graph to show all three scales by utility
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>% 
  pivot_longer(c(wLeaksRRBG, wLeaksRRCG, wLeaksRREV, 
                 wLeaksRRFG, wLeaksRRLU, wLeaksRRNG), 
               names_to = "Utility", values_to = "leakDensityBG") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races",
                       "MA_MINORITY17", "MA_INCOME17")) %>% 
  drop_na(leakDensityBG) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income"),
         Utility = recode(Utility, "wLeaksRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeaksRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeaksRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeaksRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeaksRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeaksRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, Utility),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, Utility, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensityJoinedU_tr2 <- ppLeakDensityJoinedU_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>% 
  pivot_longer(c(wLeaksRRBG, wLeaksRRCG, wLeaksRREV, 
                 wLeaksRRFG, wLeaksRRLU, wLeaksRRNG), 
               names_to = "Utility", values_to = "leakDensityTR") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityTR) %>%
  mutate(Utility = recode(Utility, "wLeaksRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeaksRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeaksRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeaksRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeaksRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeaksRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, Utility),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, Utility, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensityJoinedU_co2 <- ppLeakDensityJoinedU_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(c(wLeaksRRBG, wLeaksRRCG, wLeaksRREV, 
                 wLeaksRRFG, wLeaksRRLU, wLeaksRRNG), 
               names_to = "Utility", values_to = "leakDensityCO") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityCO) %>%
  mutate(Utility = recode(Utility, "wLeaksRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeaksRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeaksRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeaksRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeaksRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeaksRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, Utility),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, Utility, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensityJoinedU_bg2,.)

# join them all together
ppLeakDensityJoinedU_all <- list(ppLeakDensityJoinedU_tr2,
                                 ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) mutate(x, UtilityGroup = paste(Utility,Group))) %>%
  reduce(., full_join, by = "UtilityGroup") %>% 
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensityJoinedU_all2 <- list(ppLeakDensityJoinedU_tr2,
                                  ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

ppLeakDensityJoinedU_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensityJoinedU_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ Utility, drop = FALSE) +
  # facet_wrap(~ Utility, scales = "free_y", drop = FALSE) +
  labs(x = NULL, 
       y = expression(paste("Ratio of group population-weighted mean leak density (leaks/", km^2, ")", " to total population-weighted mean",sep = "")),
       title = "Relative Exposure to Unrepaired Gas Leaks in 2019 across Massachusetts by Utility", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotGraphallscalesUtilityRE.jpeg", dpi = 1000)
```
Berkshire Gas, which serves the western-most region of the state, stands out from all other utilities in Massachusetts because race and ethnicity are the least indicative of higher exposure within that utility’s service area. Indeed, Berkshire Gas is the only utility in which Whites are exposed to leak densities similar to or greater than the general population. For the Berkshire Gas territory, residents in low income environmental justice communities (MA Low Income), less educated adults, limited English speaking households, and renters, are the most disproportionately exposed groups. The ordering of those most affected groups varies slightly depending on the scale of analysis. In a few cases, the scale of analysis has a significant effect on the outcome of the analysis. In the Berkshire Gas utility territory, Blacks, Hispanics, and Housing Burdened populations may have exposures above or below that of the general population, depending on the scale of analysis. The same is true for Asians in the Liberty Utilities and Unitil/Fitchburg Gas territories. In most cases, however, there is considerable overlap or agreement between the scales of analysis.


### Statewide exposure normalized by occupied housing units

Previous research has found that distribution-level natural gas leaks occur more frequently where there are more gas service lines [@Gallagher_Down_2015; @Scott_Scott_2019]. Thus, to determine if leak frequency exceeds what would be expected given the infrastructure, leak counts should be normalized per mile of service lines. While the number or mileage of service lines in an area is available at the utility scale, this information is not available at smaller scales, such as by municipality or Census statistical units. We use the number of occupied housing units in an area as a proxy for service lines, assuming each occupied housing unit that utilizes utility-supplied natural gas requires a service line.

The geographic pattern of gas leaks per occupied housing unit differs from that of leak density (see Figures \@ref(fig:mapOHUbg) - \@ref(fig:mapOHUco)). While gas leak density is highest within the state’s most populous cities, leaks per occupied housing unit are generally highest in adjacent suburban communities, such as those to the west and south of Boston, and north of Lawrence. These differences are reflected in population-weighted relative exposures to leaks per occupied housing unit. 

```{r mapOHUbg, fig.align='center', fig.cap= "Massachusetts gas leaks per occupied housing unit by Census Block Group, 2019", cache=TRUE}
# create map of leaks per OHU by BG
blkgrpsLeaksOHU <- ma_blkgrps %>% 
  filter(ALLleaks_hu > 0)

m1b <- tm_shape(blkgrpsLeaksOHU, bbox = ma_state_sf) + 
  tm_fill(col = "ALLleaks_hu", style = "fisher",
          palette = "YlOrRd", title = "Leaks per Occupied\nHousing Unit",
          legend.format = list(digits = 2)) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_blkgrps) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.1),
            title = "Massachusetts Gas Leaks per Occupied Housing Unit by Census Block Group, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m1b, filename = "Images/pub/m_blkgrpOHU.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_blkgrpOHU.jpeg")
```

```{r mapOHUtr, fig.align='center', fig.cap= "Massachusetts gas leaks per occupied housing unit by Census Tract, 2019", cache=TRUE}
# create map of leaks per OHU by Tract
tractsLeaksOHU <- ma_tracts %>% 
  filter(ALLleaks_hu > 0)

m2b <- tm_shape(tractsLeaksOHU, bbox = ma_state_sf) + 
  tm_fill(col = "ALLleaks_hu", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Leaks per Occupied\nHousing Unit",
                                      legend.format = list(digits = 2)) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_tracts) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.10),
            title = "Massachusetts Gas Leaks per Occupied Housing Unit by Census Tract, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m2b, filename = "Images/pub/m_tractOHU.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_tractOHU.jpeg")
```

```{r mapOHUco, fig.align='center', fig.cap= "Massachusetts gas leaks per occupied housing unit by municipality, 2019", cache=TRUE}
# create map of leaks per OHU by Cosub
cosubsLeaksOHU <- ma_cosub %>% 
  filter(ALLleaks_hu > 0)

m2b <- tm_shape(cosubsLeaksOHU, bbox = ma_state_sf) + 
  tm_fill(col = "ALLleaks_hu", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Leaks per Occupied\nHousing Unit",
                                      legend.format = list(digits = 2)) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_cosub) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.10),
            title = "Massachusetts Gas Leaks per Occupied Housing Unit by Municipality, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m2b, filename = "Images/pub/m_cosubOHU.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_cosubOHU.jpeg")
```

Across all leak classes, Asians have consistently the highest relative exposures to population-weighted mean unrepaired leaks per occupied housing unit (see Figure \@ref(fig:figDotOHU)). For all leak classes combined, Whites are also slightly more exposed to higher unrepaired leaks per occupied housing unit relative to the general population.  However, the latter pattern only holds for Class 3 leaks. For the potentially more hazardous Class 1 and Class 2 leaks, Whites are the least exposed, while People of Color (especially Asians and Blacks) and limited-English speaking households are the most exposed at all scales of analysis.

```{r figDotOHU, fig.align='center', fig.cap="Relative exposures to population-weighted mean unrepaired gas leaks per occupied housing unit at Census Block Group, Census Tract, and municipality scales"}
# dot graph of relative exposure for leaks per OHU around zero line
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

# Set up for dot graph to show all three scales
ppLeakDensity_df_bg2 <- ppLeakDensity_df_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>% 
  pivot_longer(wLeaksPerHURR:wLeaksPerHURRC3, 
               names_to = "leakClass", values_to = "leakDensityBG") %>% 
  # filter(!Group %in% c("MA_ENGLISH","MA_MINORITY21","MA_INCOME21")) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income")) %>%
  mutate(leakClass = recode(leakClass, "wLeaksPerHURR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeaksPerHURRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeaksPerHURRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeaksPerHURRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, leakClass),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, leakClass, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensity_df_tr2 <- ppLeakDensity_df_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>% 
  pivot_longer(wLeaksPerHURR:wLeaksPerHURRC3, 
               names_to = "leakClass", values_to = "leakDensityTR") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wLeaksPerHURR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeaksPerHURRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeaksPerHURRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeaksPerHURRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, leakClass),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, leakClass, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensity_df_co2 <- ppLeakDensity_df_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(wLeaksPerHURR:wLeaksPerHURRC3, 
               names_to = "leakClass", values_to = "leakDensityCO") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wLeaksPerHURR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeaksPerHURRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeaksPerHURRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeaksPerHURRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, leakClass),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, leakClass, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensity_df_bg2 <- ppLeakDensity_df_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensity_df_bg2,.)

# join them all together
ppLeakDensity_df_all <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                             ppLeakDensity_df_co2) %>% 
  lapply(., function(x) mutate(x, ClassGroup = paste(leakClass,Group))) %>%
  reduce(., full_join, by = "ClassGroup") %>% 
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensity_df_all2 <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                              ppLeakDensity_df_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
ppLeakDensity_df_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensity_df_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ leakClass, scales = "free_y") +
  labs(x = NULL, 
       y = "Ratio of group population-weighted mean leaks per occupied housing unit to population mean",
       title = "Relative Exposure to Unrepaired Leaks Per Occupied Housing Unit in 2019 across Massachusetts", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotGraphOHUallscalesRE.jpeg", dpi = 1000)
```


### Utility-specific exposure normalized by occupied housing units

When broken out by utility for all leak classes combined, relative exposures are highest for non-White and low income groups in all territories except for Liberty Utilities (see Figure \@ref(fig:figDotOHUUtility)). As with the statewide pattern, however, the pattern of disparate exposure is reinforced for more hazardous leak classes. For example, Figure \@ref(fig:figDotOHUUtilityC2) shows relative exposures per occupied housing unit for Class 2 leaks by utility. In the latter case, residents in Limited English environmental justice communities within the Liberty Utilities territory have relative exposures six times that of the general population and exceed the relative exposures of all other groups. Note that Unitil/Fitchburg Gas reported no unrepaired Class 2 gas leaks in 2019.

```{r figDotOHUUtility, fig.align='center', fig.cap="Relative exposures to population-weighted mean unrepaired gas leaks per occupied housing unit by utility at Census Block Group, Census Tract, and municipality scales"}
# Set up for dot graph to show all three scales by utility
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>%
  pivot_longer(c(wLeaksPerHURRBG, wLeaksPerHURRCG, wLeaksPerHURREV, 
                 wLeaksPerHURRFG, wLeaksPerHURRLU, wLeaksPerHURRNG), 
               names_to = "Utility", values_to = "leakDensityBG") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races",
                       "MA_MINORITY17", "MA_INCOME17")) %>% 
  drop_na(leakDensityBG) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income"),
         Utility = recode(Utility, "wLeaksPerHURRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeaksPerHURRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeaksPerHURREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeaksPerHURRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeaksPerHURRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeaksPerHURRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, Utility),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, Utility, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensityJoinedU_tr2 <- ppLeakDensityJoinedU_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>% 
  pivot_longer(c(wLeaksPerHURRBG, wLeaksPerHURRCG, wLeaksPerHURREV, 
                 wLeaksPerHURRFG, wLeaksPerHURRLU, wLeaksPerHURRNG), 
               names_to = "Utility", values_to = "leakDensityTR") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityTR) %>%
  mutate(Utility = recode(Utility, "wLeaksPerHURRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeaksPerHURRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeaksPerHURREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeaksPerHURRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeaksPerHURRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeaksPerHURRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, Utility),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, Utility, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensityJoinedU_co2 <- ppLeakDensityJoinedU_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(c(wLeaksPerHURRBG, wLeaksPerHURRCG, wLeaksPerHURREV, 
                 wLeaksPerHURRFG, wLeaksPerHURRLU, wLeaksPerHURRNG), 
               names_to = "Utility", values_to = "leakDensityCO") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityCO) %>%
  mutate(Utility = recode(Utility, "wLeaksPerHURRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeaksPerHURRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeaksPerHURREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeaksPerHURRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeaksPerHURRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeaksPerHURRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, Utility),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, Utility, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensityJoinedU_bg2,.)

# join them all together
ppLeakDensityJoinedU_all <- list(ppLeakDensityJoinedU_tr2,
                                 ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) mutate(x, UtilityGroup = paste(Utility,Group))) %>%
  reduce(., full_join, by = "UtilityGroup") %>% 
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensityJoinedU_all2 <- list(ppLeakDensityJoinedU_tr2,
                                  ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

ppLeakDensityJoinedU_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensityJoinedU_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ Utility, drop = FALSE) +
  # facet_wrap(~ Utility, scales = "free_y", drop = FALSE) +
  labs(x = NULL, 
       y = "Ratio of group population-weighted mean leaks per occupied housing unit to population mean",
       title = "Relative Exposure to Unrepaired Leaks Per Occupied Housing Unit in 2019 across Massachusetts by Utility", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotGraphOHUallscalesUtilityRE.jpeg", dpi = 1000)
```

```{r figDotOHUUtilityC2, fig.align='center', fig.cap="Relative exposures to population-weighted mean unrepaired class 2 gas leaks per occupied housing unit by utility at Census Block Group, Census Tract, and municipality scales"}
# Set up for dot graph to show all three scales by utility
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>%
  pivot_longer(c(wLeaksPerHURRBGC2, wLeaksPerHURRCGC2, wLeaksPerHURREVC2, 
                 wLeaksPerHURRFGC2, wLeaksPerHURRLUC2, wLeaksPerHURRNGC2), 
               names_to = "Utility", values_to = "leakDensityBG") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races",
                       "MA_MINORITY17", "MA_INCOME17")) %>% 
  drop_na(leakDensityBG) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income"),
         Utility = recode(Utility, "wLeaksPerHURRBGC2" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRCGC2" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURREVC2" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRLUC2" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRNGC2" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRFGC2" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg" & Class == "2")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, Utility),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, Utility, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensityJoinedU_tr2 <- ppLeakDensityJoinedU_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>% 
  pivot_longer(c(wLeaksPerHURRBGC2, wLeaksPerHURRCGC2, wLeaksPerHURREVC2, 
                 wLeaksPerHURRFGC2, wLeaksPerHURRLUC2, wLeaksPerHURRNGC2), 
               names_to = "Utility", values_to = "leakDensityTR") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityTR) %>%
  mutate(Utility = recode(Utility, "wLeaksPerHURRBGC2" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRCGC2" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURREVC2" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRLUC2" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRNGC2" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRFGC2" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg" & Class == "2")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, Utility),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, Utility, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensityJoinedU_co2 <- ppLeakDensityJoinedU_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(c(wLeaksPerHURRBGC2, wLeaksPerHURRCGC2, wLeaksPerHURREVC2, 
                 wLeaksPerHURRFGC2, wLeaksPerHURRLUC2, wLeaksPerHURRNGC2), 
               names_to = "Utility", values_to = "leakDensityCO") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityCO) %>%
  mutate(Utility = recode(Utility, "wLeaksPerHURRBGC2" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRCGC2" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURREVC2" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRLUC2" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRNGC2" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid" & Class == "2")),big.mark = ","),")"),
                          "wLeaksPerHURRFGC2" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg" & Class == "2")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, Utility),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, Utility, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensityJoinedU_bg2,.)

# join them all together
ppLeakDensityJoinedU_all <- list(ppLeakDensityJoinedU_tr2,
                                 ppLeakDensityJoinedU_bg2,
                                 ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) mutate(x, UtilityGroup = paste(Utility,Group))) %>%
  reduce(., full_join, by = "UtilityGroup") %>% 
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensityJoinedU_all2 <- list(ppLeakDensityJoinedU_tr2,
                                  ppLeakDensityJoinedU_bg2,
                                  ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

ppLeakDensityJoinedU_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensityJoinedU_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ Utility, drop = FALSE) +
  # facet_wrap(~ Utility, scales = "free_y", drop = FALSE) +
  labs(x = NULL, 
       y = "Ratio of group population-weighted mean leaks per occupied housing unit to population mean",
       title = "Relative Exposure to Unrepaired Class 2 Leaks Per Occupied Housing Unit in 2019 across Massachusetts by Utility", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotGraphOHUallscalesC2UtilityRE.jpeg", dpi = 1000)
```


## Response – Leak Age and Repair
Another way to evaluate the equity of natural gas leaks is to consider how responses to those leaks vary. This section considers two measures of leak response: 

* the age of leaks when they were repaired in 2019, and 
* the age of leaks that remained unrepaired at the end of the 2019 reporting year. 


### Age of repaired leaks

In 2019, repairs were reported for `r format(nrow(repaired2019),big.mark = ",")` leaks (`r total_Sumdf %>% filter(Class == "Total") %>% select(PctRepaired) %>% pull() %>% round(., 1)`% of all reported leaks) across the state (see Table \@ref(tab:tabREP) and Figure \@ref(fig:histREP)). More than half of all leaks repaired were Class 1 leaks (`r format(repaired2019 %>% filter(Class == "1") %>% nrow(),big.mark = ",")` or `r round(total_Sumdf[1,"Repaired"]/total_Sumdf[4,"Repaired"]*100,0)`%), which is consistent with state policy which requires leaks of greatest potential hazard to be repaired immediately. Class 2 leaks were `r format(repaired2019 %>% filter(Class == "2") %>% nrow(),big.mark = ",")` leaks (`r round(total_Sumdf[2,"Repaired"]/total_Sumdf[4,"Repaired"]*100,0)`% of leaks repaired), followed distantly by Class 3 leaks, which constituted only `r format(repaired2019 %>% filter(Class == "3") %>% nrow(),big.mark = ",")` leaks (`r round(total_Sumdf[3,"Repaired"]/total_Sumdf[4,"Repaired"]*100,0)`% of all repaired leaks in 2019). Age for repaired leaks was calculated as the duration of time, in days, from the date the leak was reported to the date of its reported repair. Ages of repaired leaks varied dramatically both within and between leak classes and across utilities.  
\pagebreak
```{r tabREP, fig.align='center'}
# time to repair by class
repaired_Timedf <- repaired2019 %>% 
  as.data.frame() %>% 
  select(Class,DaysToRepair) %>% 
  drop_na(Class) %>% 
  group_by(Class) %>% 
  summarize(across(.cols = DaysToRepair, summary_stats, 
                   .names = "{.fn}")) %>% 
  arrange(Class)

# add row for all repaired
repaired_Timedf <- repaired2019 %>% 
  as.data.frame() %>% 
  select(DaysToRepair) %>% 
  summarize(across(.cols = DaysToRepair, summary_stats, 
                   .names = "{.fn}")) %>% 
  mutate(Class = "All", .before = Min) %>% 
  rbind(repaired_Timedf, .)

repaired_Timedf %>% kable(., longtable = T, booktabs = T,
                          format.args = list(big.mark = ','),
                          caption = "Repaired Leak Age (days), 2019", 
                          align = "r", 
                          digits = 1) %>% 
  # add_header_above(., c(" ", "Unrepaired Leaks" = 4, "Repaired Leaks" = 4, 
  #                       "Total Leaks" = 4)) %>% 
  # row_spec(4, bold = T) %>% 
  kable_styling(latex_options = c("repeat_header"))
```

```{r histREP, fig.align='center', fig.cap="Histogram of ages of leaks repaired in 2019 by leak class"}
# create group medians for histogram
mu1 <- repaired2019 %>%
  as.data.frame() %>%
  drop_na(Class) %>%
  group_by(Class) %>%
  summarize(grp.med = median(DaysToRepair, na.rm = T))

repaired2019 %>% 
  as.data.frame() %>% 
  select(Class,DaysToRepair) %>% 
  drop_na(Class) %>% 
  ggplot(aes(x = DaysToRepair, fill = Class, color = Class)) +
  geom_histogram(alpha = 0.6, position = "identity") +
  geom_vline(data = mu1, aes(xintercept = grp.med, color = Class), 
             linetype = "dashed") +
  geom_hline(yintercept = 1, color = "gray88") +
  scale_color_manual(values = c("#fc8d62","#66c2a5","#8da0cb"),
                     guide = guide_legend(override.aes = list(fill = c("#fc8d62","#66c2a5","#8da0cb")))) +
  scale_fill_manual(values = c("#fc8d62","#66c2a5","#8da0cb"),
                    guide = FALSE) +
  theme_minimal() +
  theme(plot.caption = element_text(hjust = 0)) +
  labs(title = "Repaired Leak Age by Leak Class in 2019",
       x = "Days to Repair",
       y = "Number of Repaired Leaks", fill = "",
       caption = "Column heights indicate number of leaks for a given time range. Colors correspond to leak class.\nDashed lines indicate median times by leak class. Rug plot at base shows individual leak ages.") +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) + geom_rug()

ggsave("Images/pub/histRepairTime.jpeg", dpi = 1000)
```

There is no clear geographic pattern or clustering of areas with older repaired leaks (see Figures \@ref(fig:mapTimeBG) - \@ref(fig:mapTimeCO)). Indeed, communities with the oldest repaired leaks are scattered throughout the state. Notably, these oldest repaired leaks are not concentrated within major urban areas but are instead in outlying suburbs.

```{r mapTimeBG, fig.align='center', fig.cap= "Average age of leaks repaired in 2019 across Massachusetts by Census Block Group", cache=TRUE}
# create a map of average repair times by block group
blkgrpsLeaksTime <- ma_blkgrps %>% 
  filter(repaired2019total > 0)

m4 <- tm_shape(blkgrpsLeaksTime, bbox = ma_state_sf) + 
  tm_fill(col = "DaysToRepairAvg", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Mean Leak\nAge (days)",
                                      legend.format = list(digits = 0),
                                      colorNA = NULL) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_blkgrps) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.1),
            title = "Massachusetts Mean Age of Repaired Leaks by Census Block Group and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)
# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m4, filename = "Images/pub/m_blkgrpTime.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_blkgrpTime.jpeg")
```

```{r mapTimeTR, fig.align='center', fig.cap= "Average age of leaks repaired in 2019 across Massachusetts by Census Tract", cache=TRUE}
# create a map of average repair times by tract
tractsLeaksTime <- ma_tracts %>% 
  filter(repaired2019total > 0)

m4 <- tm_shape(tractsLeaksTime, bbox = ma_state_sf) + 
  tm_fill(col = "DaysToRepairAvg", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Mean Leak\nAge (days)",
                                      legend.format = list(digits = 0),
                                      colorNA = NULL) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_tracts) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.10),
            title = "Massachusetts Mean Age of Repaired Leaks by Census Tract and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m4, filename = "Images/pub/m_tractTime.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_tractTime.jpeg")
```

```{r mapTimeCO, fig.align='center', fig.cap= "Average age of leaks repaired in 2019 across Massachusetts by municipality", cache=TRUE}
# create a map of average repair times by cosub
cosubLeaksTime <- ma_cosub %>% 
  filter(repaired2019total > 0)

m4 <- tm_shape(cosubLeaksTime, bbox = ma_state_sf) + 
  tm_fill(col = "DaysToRepairAvg", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Mean Leak\nAge (days)",
                                      legend.format = list(digits = 0),
                                      colorNA = NULL) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_cosub) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.10),
            title = "Massachusetts Mean Age of Repaired Leaks by Municipality and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m4, filename = "Images/pub/m_cosubTime.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_cosubTime.jpeg")
```

For all leak classes combined, People of Color (especially Blacks and Asians) and renters experience the longest repair times (i.e., higher relative exposures to population-weighted mean age of repaired leaks compared to the general population) (see Figure \@ref(fig:figDotTime)). This pattern holds true for Class 2 and Class 3 leaks. For Class 1 leaks, residents of state-designated low income environmental justice communities have the greatest relative exposure. Whites also experience consistently greater population-weighted age for repaired Class 1 leaks.  Note however that  the differences in age between Class 1 leaks average only a few days, while differences in average age for Class 2 and 3 leaks are more than a year (see Table \@ref(tab:tabREP)). 

```{r figDotTime, fig.align='center', fig.cap="Relative exposures to population-weighted mean age of leaks repaired in 2019 at Census Block Group, Census Tract, and municipality scales"}
# dot graph of relative exposure for unrepaired leak density around zero line
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

# Set up for dot graph to show all three scales
ppLeakDensity_df_bg2 <- ppLeakDensity_df_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>%
  pivot_longer(wDaysToRepairAvgRR:wDaysToRepairAvgRRC3, 
               names_to = "leakClass", values_to = "leakDensityBG") %>% 
  # filter(!Group %in% c("MA_ENGLISH","MA_MINORITY21","MA_INCOME21")) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income")) %>%
  mutate(leakClass = recode(leakClass, "wDaysToRepairAvgRR" = paste0("All Leaks"," (n = ",formatC(nrow(repaired2019),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "1")),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "2")),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, leakClass),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, leakClass, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensity_df_tr2 <- ppLeakDensity_df_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>%
  pivot_longer(wDaysToRepairAvgRR:wDaysToRepairAvgRRC3, 
               names_to = "leakClass", values_to = "leakDensityTR") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wDaysToRepairAvgRR" = paste0("All Leaks"," (n = ",formatC(nrow(repaired2019),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "1")),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "2")),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, leakClass),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, leakClass, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensity_df_co2 <- ppLeakDensity_df_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>% 
  pivot_longer(wDaysToRepairAvgRR:wDaysToRepairAvgRRC3, 
               names_to = "leakClass", values_to = "leakDensityCO") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wDaysToRepairAvgRR" = paste0("All Leaks"," (n = ",formatC(nrow(repaired2019),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "1")),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "2")),big.mark = ","),")"),
                            "wDaysToRepairAvgRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(repaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, leakClass),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, leakClass, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensity_df_bg2 <- ppLeakDensity_df_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensity_df_bg2,.)

# join them all together
ppLeakDensity_df_all <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                             ppLeakDensity_df_co2) %>% 
  lapply(., function(x) mutate(x, ClassGroup = paste(leakClass,Group))) %>%
  reduce(., full_join, by = "ClassGroup") %>%  
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensity_df_all2 <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                              ppLeakDensity_df_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
ppLeakDensity_df_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensity_df_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ leakClass, scales = "free_y") +
  labs(x = NULL, 
       y = "Ratio of group population-weighted mean leak repair time to population mean",
       title = "Relative Exposure to Mean Age of Leaks Repaired in 2019 across Massachusetts", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotGraphTimeallscalesRE.jpeg", dpi = 1000)
```


### Utility-specific ages of repaired leaks

Median ages of repaired leaks vary only slightly by utility, but there are significant differences in terms of outliers (see Figure \@ref(fig:boxTimeU); graphs of utility-specific leak ages by leak class are available in supplementary materials). Liberty Utilities and Berkshire Gas exhibited the highest median age of repaired leaks, but National Grid showed the highest frequency and magnitudes of outliers. 

```{r boxTimeU, fig.align='center', fig.cap="Boxplot of ages of repaired gas leaks by utility in 2019"}
# # ordered boxplot of time to repair by utility
# repaired2019 %>% 
#   as.data.frame() %>% 
#   mutate(Utility = recode(Utility, "National Grid - Boston Gas" = "National Grid",
#                           "National Grid - Colonial Gas" = "National Grid",
#                           "Fitchburg" = "Unitil/Fitchburg Gas")) %>% 
#   ggplot(aes(x = reorder(Utility, DaysToRepair, FUN = median, na.rm = TRUE), 
#              y = DaysToRepair, fill = Utility)) +
#   # geom_violin(width=1.4) +
#   # geom_boxplot(width=0.1, color="grey", alpha=0.2) + 
#   geom_boxplot() + 
#   coord_flip() +
#   theme_minimal() +
#   theme(
#     legend.position="none") +
#   ggtitle("Ages (days) of Repaired Leaks by Utility in 2019") +
#   xlab("") + ylab("Days to repair through Dec 31, 2019") + 
#   scale_y_continuous(labels = function(x) format(x, big.mark = ","))
#   # facet_wrap("Class", scales = "free")
# 
# ggsave(filename = "Images/pub/boxplotLeaksByUtilityTime.jpeg")


# create a faceted boxplot for utility leak repair ages by leak class
# create simplified dfs to rbind 
boxplotUtime <- repaired2019 %>% 
  as.data.frame() %>% 
  transmute(Class = recode(Class, 
                           "1" = paste0("Class 1 Leaks",
                                        " (n = ", 
                                        formatC(nrow(filter(., Class == "1")),big.mark = ","),")"),
                           "2" = paste0("Class 2 Leaks",
                                        " (n = ", 
                                        formatC(nrow(filter(., Class == "2")),big.mark = ","),")"),
                           "3" = paste0("Class 3 Leaks",
                                        " (n = ", 
                                        formatC(nrow(filter(., Class == "3")),big.mark = ","),")")),
            Utility = recode(Utility, 
                          "National Grid - Boston Gas" = "National Grid",
                          "National Grid - Colonial Gas" = "National Grid",
                          "Fitchburg" = "Unitil/Fitchburg Gas"),
            DaysToRepair = DaysToRepair)

# create "all leak class", assemble df and generate faceted graph
boxplotUtime %>% 
  mutate(Class = paste0("All Leaks", " (n = ", 
                                        formatC(nrow(.),
                                                big.mark = ","),")")) %>% 
  rbind(., boxplotUtime) %>% 
  ggplot(aes(x = reorder(Utility, DaysToRepair, FUN = median, 
                         na.rm = TRUE), 
             y = DaysToRepair, fill = Utility)) +
  geom_boxplot(outlier.size = 0.9) + 
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position="none") +
  ggtitle("Ages (days) of Repaired Leaks by Utility in 2019") +
  xlab("") + ylab("Days to repair through Dec 31, 2019") + 
  scale_y_continuous(labels = function(x) format(x, big.mark = ",")) +
  facet_wrap("Class")

ggsave(filename = "Images/pub/boxplotLeaksByUtilityTime.jpeg", dpi = 1000)
```

People of Color (especially Blacks and Asians), followed by lower income groups, are most exposed to the oldest repaired leaks in the state’s three largest utilities – National Grid, Eversource Energy, and Columbia Gas (see Figure \@ref(fig:figDotTimeU)). Interestingly, Whites are also consistently exposed to repaired leak ages slightly above the general population in four out of the six utility territories, albeit primarily in the smaller utilities and never as the most exposed group. In utility territories where People of Color are less exposed, it is instead groups with less education (Berkshire Gas), people who are older (Berkshire Gas, Liberty Utilities), and residents of lower income environmental justice communities (Unitil/Fitchburg Gas) who are exposed to the oldest repaired leak ages. 

```{r figDotTimeU, fig.align='center', fig.cap="Relative exposures to population-weighted mean age of leaks repaired in 2019 by utility at Census Block Group, Census Tract, and municipality scales"}
# Set up for dot graph to show all three scales by utility
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>% 
  pivot_longer(c(wDaysToRepairAvgRRBG, wDaysToRepairAvgRRCG, wDaysToRepairAvgRREV, 
                 wDaysToRepairAvgRRFG, wDaysToRepairAvgRRLU, wDaysToRepairAvgRRNG), 
               names_to = "Utility", values_to = "leakDensityBG") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races",
                       "MA_MINORITY17", "MA_INCOME17")) %>% 
  drop_na(leakDensityBG) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income"),
         Utility = recode(Utility, "wDaysToRepairAvgRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wDaysToRepairAvgRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(repaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, Utility),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, Utility, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensityJoinedU_tr2 <- ppLeakDensityJoinedU_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>%
  pivot_longer(c(wDaysToRepairAvgRRBG, wDaysToRepairAvgRRCG, wDaysToRepairAvgRREV, 
                 wDaysToRepairAvgRRFG, wDaysToRepairAvgRRLU, wDaysToRepairAvgRRNG), 
               names_to = "Utility", values_to = "leakDensityTR") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityTR) %>%
  mutate(Utility = recode(Utility, "wDaysToRepairAvgRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wDaysToRepairAvgRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(repaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, Utility),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, Utility, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensityJoinedU_co2 <- ppLeakDensityJoinedU_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(c(wDaysToRepairAvgRRBG, wDaysToRepairAvgRRCG, wDaysToRepairAvgRREV, 
                 wDaysToRepairAvgRRFG, wDaysToRepairAvgRRLU, wDaysToRepairAvgRRNG), 
               names_to = "Utility", values_to = "leakDensityCO") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityCO) %>%
  mutate(Utility = recode(Utility, "wDaysToRepairAvgRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wDaysToRepairAvgRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(repaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wDaysToRepairAvgRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(repaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, Utility),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, Utility, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensityJoinedU_bg2,.)

# join them all together
ppLeakDensityJoinedU_all <- list(ppLeakDensityJoinedU_tr2, ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) mutate(x, UtilityGroup = paste(Utility,Group))) %>%
  reduce(., full_join, by = "UtilityGroup") %>%  
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensityJoinedU_all2 <- list(ppLeakDensityJoinedU_tr2,
                                  ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

ppLeakDensityJoinedU_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensityJoinedU_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ Utility, drop = FALSE) +
  # facet_wrap(~ Utility, scales = "free_y", drop = FALSE) +
  labs(x = NULL, 
       y = "Ratio of group population-weighted mean leak repair time(days) to general population-weighted mean",
       title = "Relative Exposure to Average Age of Leaks Repaired in 2019 across Massachusetts by Utility", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotTimeUtilityRE.jpeg", dpi = 1000)
```


### Age of unrepaired leaks

```{r tabAgeCalc,include=FALSE}
# unrepaired age by class
unrepaired_Timedf <- unrepaired2019 %>% 
  as.data.frame() %>% 
  select(Class,LeakAgeDays) %>% 
  drop_na(Class) %>% 
  group_by(Class) %>% 
  summarize(across(.cols = LeakAgeDays, summary_stats, 
                                       .names = "{.fn}")) %>% 
  arrange(Class)

# add row for all unrepaired
unrepaired_Timedf <- unrepaired2019 %>% 
  as.data.frame() %>% 
  select(LeakAgeDays) %>% 
  summarize(across(.cols = LeakAgeDays, summary_stats, 
                   .names = "{.fn}")) %>% 
  mutate(Class = "All", .before = Min) %>% 
  rbind(unrepaired_Timedf, .)
```

At the end of calendar year 2019, there were `r format(nrow(unrepaired2019),big.mark = ",")` unrepaired natural gas leaks across the state, representing `r total_Sumdf %>% filter(Class == "Total") %>% select(PctUnrepaired) %>% pull() %>% round(., 1)`% of all leaks reported in that year (see Table \@ref(tab:tabAge) and Figure \@ref(fig:histAge)). Most of the unrepaired leaks were Class 3 leaks (`r round(total_Sumdf[3,"Unrepaired"]/total_Sumdf[4,"Unrepaired"]*100,0)`%), the least hazardous leak classification with the most lenient repair protocols. Class 2 leaks were `r round(total_Sumdf[2,"Unrepaired"]/total_Sumdf[4,"Unrepaired"]*100,1)`% of unrepaired leaks, followed distantly by Class 1 leaks, which constituted only `r round(total_Sumdf[1,"Unrepaired"]/total_Sumdf[4,"Unrepaired"]*100,2)`% (`r unrepaired2019 %>% filter(Class == "1") %>% nrow()`) of all unrepaired leaks at the end of 2019. The age of these unrepaired leaks was calculated as the duration of time from the date the leak was reported to the end of the calendar year – December 31, 2019. 

Leak ages varied dramatically both within and between leak classes and across utilities. The median age of Class 1 leaks was `r unrepaired_Timedf["1","Med"] %>% pull()` days, `r unrepaired_Timedf["2","Med"] %>% pull()` days (`r round(unrepaired_Timedf["2","Med"]/30,1)` months) for Class 2 leaks, and `r format(unrepaired_Timedf["3","Med"] %>% pull(), big.mark = ",")` days (`r round(unrepaired_Timedf["3","Med"]/365,1)` years) for Class 3 leaks.

```{r tabAge, fig.align='center'}
# unrepaired age by class
unrepaired_Timedf <- unrepaired2019 %>% 
  as.data.frame() %>% 
  select(Class,LeakAgeDays) %>% 
  drop_na(Class) %>% 
  group_by(Class) %>% 
  summarize(across(.cols = LeakAgeDays, summary_stats, 
                                       .names = "{.fn}")) %>% 
  arrange(Class)

# add row for all unrepaired
unrepaired_Timedf <- unrepaired2019 %>% 
  as.data.frame() %>% 
  select(LeakAgeDays) %>% 
  summarize(across(.cols = LeakAgeDays, summary_stats, 
                   .names = "{.fn}")) %>% 
  mutate(Class = "All", .before = Min) %>% 
  rbind(unrepaired_Timedf, .)

unrepaired_Timedf %>% kable(., longtable = T, booktabs = T,
                          format.args = list(big.mark = ','),
                          caption = "Unrepaired gas leak age (days), 2019", 
                          align = "r", 
                          digits = 1) %>% 
  # add_header_above(., c(" ", "Unrepaired Leaks" = 4, "Repaired Leaks" = 4, 
  #                       "Total Leaks" = 4)) %>% 
  # row_spec(4, bold = T) %>% 
  kable_styling(latex_options = c("repeat_header"))
```

```{r histAge, fig.align='center', fig.cap="Histogram of unrepaired leak age by leak class"}
# create group medians for histogram
mu <- unrepaired2019 %>% 
  as.data.frame() %>% 
  drop_na(Class) %>% 
  group_by(Class) %>% 
  summarize(grp.med = median(LeakAgeDays, na.rm = T))

unrepaired2019 %>% 
  as.data.frame() %>% 
  select(Class,LeakAgeDays) %>% 
  mutate(Class = factor(Class, levels = c("3","2","1"), ordered = T)) %>% 
  drop_na(Class) %>% 
  ggplot(aes(x = LeakAgeDays, fill = Class, color = Class)) +
  geom_histogram(alpha = 0.6, position = "identity") +
  geom_vline(data = mu, aes(xintercept = grp.med, color = Class), 
             linetype = "dashed") +
  geom_hline(yintercept = 0, color = "grey88") +
  scale_color_manual(values = c("#8da0cb","#66c2a5","#fc8d62"),
                     guide = guide_legend(reverse = TRUE,
                                          override.aes = list(fill = c("#fc8d62","#66c2a5","#8da0cb")))) +
  scale_fill_manual(values = c("#8da0cb","#66c2a5","#fc8d62"),
                    guide = FALSE) +
  theme_minimal() +
  theme(plot.caption = element_text(hjust = 0)) +
  labs(title = "Unrepaired Leak Age by Leak Class in 2019",
       x = "Age (days)",
       y = "Number of Unrepaired Leaks", fill = "",
       caption = "Column heights indicate number of leaks for a given age range. Colors correspond to leak class.\nDashed lines indicate median age by leak class. Rug plot at base shows individual leak ages.") +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) + geom_rug()

ggsave("Images/pub/histAge.jpeg", dpi = 1000)
```

Geographically, areas with the oldest unrepaired leaks radiate outward from Boston and inner core communities, primarily along a southeast to northwest axis (see Figures \@ref(fig:magBGage) - \@ref(fig:magCOage)). Areas with the oldest unrepaired leaks are concentrated in the eastern half of the state, with some outliers in the southeast and south central regions. 

```{r magBGage, fig.align='center', fig.cap= "Average age of unrepaired leaks in 2019 across Massachusetts by Census Block Group", cache=TRUE}
# create a map of average leak age by block group
blkgrpsLeaksAge <- ma_blkgrps %>% 
  filter(unrepaired2019total > 0)

m5 <- tm_shape(blkgrpsLeaksAge, bbox = ma_state_sf) + 
  tm_fill(col = "LeakAgeDaysAvg", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Mean Leak\nAge (days)",
                                      legend.format = list(digits = 0),
                                      colorNA = NULL) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_blkgrps) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.10),
            title = "Massachusetts Mean Age of Unrepaired Leaks by Census Block Group and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m5, filename = "Images/pub/m_blkgrpAge.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_blkgrpAge.jpeg")
```

```{r magTRage, fig.align='center', fig.cap= "Average age of unrepaired leaks in 2019 across Massachusetts by Census Tract", cache=TRUE}
# create a map of average leak age by block group
tractsLeaksAge <- ma_tracts %>% 
  filter(unrepaired2019total > 0)

m5 <- tm_shape(tractsLeaksAge, bbox = ma_state_sf) + 
  tm_fill(col = "LeakAgeDaysAvg", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Mean Leak\nAge (days)",
                                      legend.format = list(digits = 0),
                                      colorNA = NULL) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_tracts) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.10),
            title = "Massachusetts Mean Age of Unrepaired Leaks by Census Tract and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m5, filename = "Images/pub/m_tractAge.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_tractAge.jpeg")
```

```{r magCOage, fig.align='center', fig.cap= "Average age of unrepaired leaks in 2019 across Massachusetts by municipality", cache=TRUE}
# create a map of average leak age by municipality
cosubLeaksAge <- ma_cosub %>% 
  filter(unrepaired2019total > 0)

m5 <- tm_shape(cosubLeaksAge, bbox = ma_state_sf) + 
  tm_fill(col = "LeakAgeDaysAvg", style = "fisher",
                                      palette = "YlOrRd", 
                                      title = "Mean Leak\nAge (days)",
                                      legend.format = list(digits = 0),
                                      colorNA = NULL) +
  tm_shape(ng_nogas_muni) + 
  tm_fill(col = "GAS_LABEL", palette = c("ivory2","gray88"), title = "", 
          legend.show = FALSE) + 
  tm_shape(ma_cosub) + tm_borders(lwd = 0.01, alpha = 0.8) +
  tm_shape(ng_service_areas2) + tm_borders(lwd = 0.8, col = "grey68", alpha = 0.7) +
  tm_shape(ng_service_labels) + 
  tm_text("ID", size = 0.6, col = "gray") +
  tm_shape(ng_service_labels2) + 
  tm_text("ID", size = 0.4, col = "gray") +
  tm_shape(ng_service_labels3) + 
  tm_text("ID", size = 0.6, col = "gray", xmod = 2.7, ymod = 0.7) +
  # tm_shape(EJlayer) + tm_borders(col = "blue", alpha = 0.7) +
  tm_shape(ma_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(newton) + tm_dots() +
  tm_text("NAME", size = 0.4, col = "black",
          xmod = -0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0,25,50), position = c("left","BOTTOM")) +
  tm_layout(legend.position = c(0.050,0.10),
            title = "Massachusetts Mean Age of Unrepaired Leaks by Municipality and Utility Territories, 2019",
            inner.margins = c(0.02,0.02,0.09,0.02)) +
  tm_credits("BG = Berkshire Gas\nBGC = Blackstone Gas Co.\nCG = Columbia Gas\nEV = Eversource Energy\nLU = Liberty Utilities\nNG = National Grid\nUN = Unitil/Fitchburg Gas", position = c(.3,.11), col = "gray47", size = 0.7)

# create insets with dummy layers to create separate legend items
d_inset <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "ivory2", labels = "Municipal Utility", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.06), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

d_inset2 <- tm_shape(ng_nogas_muni) + tm_fill(alpha = 0) +
  tm_add_legend(type = "fill", col = "gray88", labels = "No Natural Gas Service", 
                border.col = "white") +
  tm_layout(legend.position = c(0.33,0.025), bg.color = NA, frame = FALSE, 
            legend.text.color = "gray47")

# create viewports to specify dimensions of insets
vp_inset <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)
vp_inset2 <- grid::viewport(x = 0.5, y = 0.5, width = 3, height = 1)

tmap_save(m5, filename = "Images/pub/m_cosubAge.jpeg", 
          dpi = 1000, height = 4, width = 8, units = "in",
          insets_tm = list(d_inset,d_inset2),
          insets_vp = list(vp_inset,vp_inset2))

knitr::include_graphics("Images/pub/m_cosubAge.jpeg")
```

For all leak classes combined, People of Color (especially Asians and Blacks), limited English speaking households, and renters live in areas where the population-weighted mean age of these unrepaired leaks is 6% to 19% higher than the general population (see Figure \@ref(fig:figDotAge)). For Class 3 leaks, which constitute the majority of unrepaired leaks (`r round(total_Sumdf[3,"Unrepaired"]/total_Sumdf[4,"Unrepaired"]*100,0)`%), residents in limited English speaking environmental justice communities are exposed to leak ages that are `r round((ppLeakDensity_df_bg %>% filter(Group == "MA_ENGLISH") %>% select(wLeakAgeDaysAvgRRC3) %>% pull()-1)*100,0)`% higher than the general population, while Asians are exposed to leak ages that are `r round((ppLeakDensity_df_bg %>% filter(Group == "Asian") %>% select(wLeakAgeDaysAvgRRC3) %>% pull()-1)*100,0)`% (Block Group level) to `r round((ppLeakDensity_df_tr %>% filter(Group == "Asian") %>% select(wLeakAgeDaysAvgRRC3) %>% pull()-1)*100,0)`% (Tract level) higher. The differences are greater for more hazardous Class 2 leaks, where Asians are exposed to unrepaired Class 2 leak ages that are `r round((ppLeakDensity_df_bg %>% filter(Group == "Asian") %>% select(wLeakAgeDaysAvgRRC2) %>% pull()-1)*100,0)`% (Block Group level) to `r round((ppLeakDensity_df_co %>% filter(Group == "Asian") %>% select(wLeakAgeDaysAvgRRC2) %>% pull()-1)*100,0)`% (Municipality level) higher. Class 1 leaks, which constitute a small (0.25%) and generally short-lived fraction of unrepaired leaks, are consistently older for lower income, disabled, and White residents.

```{r figDotAge, fig.align='center', fig.cap="Relative exposure to population-weighted mean age of unrepaired gas leaks in 2019 at Census Block Group, Census Tract, and municipality scales"}
# dot graph of relative exposure for unrepaired leak age around zero line
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

# Set up for dot graph to show all three scales
ppLeakDensity_df_bg2 <- ppLeakDensity_df_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>%
  pivot_longer(wLeakAgeDaysAvgRR:wLeakAgeDaysAvgRRC3, 
               names_to = "leakClass", values_to = "leakDensityBG") %>% 
  # filter(!Group %in% c("MA_ENGLISH","MA_MINORITY21","MA_INCOME21")) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income")) %>%
  mutate(leakClass = recode(leakClass, "wLeakAgeDaysAvgRR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, leakClass),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, leakClass, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensity_df_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensity_df_tr2 <- ppLeakDensity_df_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>%
  pivot_longer(wLeakAgeDaysAvgRR:wLeakAgeDaysAvgRRC3, 
               names_to = "leakClass", values_to = "leakDensityTR") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wLeakAgeDaysAvgRR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, leakClass),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, leakClass, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensity_df_co2 <- ppLeakDensity_df_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(wLeakAgeDaysAvgRR:wLeakAgeDaysAvgRRC3, 
               names_to = "leakClass", values_to = "leakDensityCO") %>% 
  # filter(!Group %in% c("Disabled Adults", "Housing Burdened")) %>%
  mutate(leakClass = recode(leakClass, "wLeakAgeDaysAvgRR" = paste0("All Leaks"," (n = ",formatC(nrow(unrepaired2019),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC1" = paste0("Class 1 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "1")),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC2" = paste0("Class 2 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "2")),big.mark = ","),")"),
                            "wLeakAgeDaysAvgRRC3" = paste0("Class 3 Leaks"," (n = ",formatC(nrow(filter(unrepaired2019,Class == "3")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, leakClass),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, leakClass, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensity_df_bg2 <- ppLeakDensity_df_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensity_df_bg2,.)

# join them all together
ppLeakDensity_df_all <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                             ppLeakDensity_df_co2) %>% 
  lapply(., function(x) mutate(x, ClassGroup = paste(leakClass,Group))) %>%
  reduce(., full_join, by = "ClassGroup") %>%  
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensity_df_all2 <- list(ppLeakDensity_df_tr2, ppLeakDensity_df_bg2,
                              ppLeakDensity_df_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
ppLeakDensity_df_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensity_df_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ leakClass, scales = "free_y") +
  labs(x = NULL, 
       y = "Ratio of group population-weighted mean leak age to population mean",
       title = "Relative Exposure to Mean Age of Unrepaired Leaks in 2019 across Massachusetts", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotGraphAgeallscalesRE.jpeg", dpi = 1000)
```


### Utility-specific ages of unrepaired leaks

Ages of unrepaired leaks varied significantly by utility (see Figure \@ref(fig:boxAgeU)). Liberty Utilities exhibited the highest median age of unrepaired leaks, while Unitil/Fitchburg Gas had the lowest. National Grid, the largest utility in the state, shows the greatest range in unrepaired leak ages. National Grid, Liberty Utilities, and Columbia Gas show a significant number of outliers for leak age. Across those three utilities, `r format(unrepaired2019 %>% filter(LeakAgeDays>(15*365)) %>% nrow(),big.mark = ",")` unrepaired leaks were over 15 years old (most of which belonged to National Grid).

When broken out by utility, non-White groups and limited English speaking households are consistently more exposed to older unrepaired leaks in four out of six utility territories, including the three largest (see Figure \@ref(fig:figDotAgeU); graphs of utility-specific leak ages by leak class are available in supplementary materials). Lower income groups and renters are also consistently exposed to older unrepaired leak ages than the general population in the territories of National Grid, Eversource, and Liberty Utilities. Berkshire Gas is the only territory where Whites are consistently the most exposed to older leak ages than any other group. Unitil/Fitchburg Gas stands out because of the wide variation in apparent exposure depending on the scale of analysis; no group is consistently (across all scales) more or less exposed. The latter result likely reflects the very small number of unrepaired leaks for that utility (`r unrepaired2019 %>% filter(Utility == "Fitchburg") %>% nrow()`) in 2019. 

```{r boxAgeU, fig.align='center', fig.cap="Boxplot of unrepaired gas leak age by utility in 2019"}
# # ordered boxplot of age of unrepaired leak ages by utility
# unrepaired2019 %>% 
#   as.data.frame() %>% 
#   mutate(Utility = recode(Utility, "National Grid - Boston Gas" = "National Grid",
#                           "National Grid - Colonial Gas" = "National Grid",
#                           "Fitchburg" = "Unitil/Fitchburg")) %>% 
#   ggplot(aes(x = reorder(Utility, LeakAgeDays, FUN = median, na.rm = TRUE), 
#              y = LeakAgeDays, fill = Utility)) +
#   # geom_violin(width=1.4) +
#   # geom_boxplot(width=0.1, color="grey", alpha=0.2) + 
#   geom_boxplot() + 
#   coord_flip() +
#   theme_minimal() +
#   theme(
#     legend.position="none") +
#   ggtitle("Age (days) of unrepaired leaks by utility") +
#   xlab("") + ylab("Days since leak reported to end of 2019") +
#   scale_y_continuous(labels = function(x) format(x, big.mark = ","))
#   # facet_wrap("Class", scales = "free")

# create a faceted boxplot for utility leak repair ages by leak class
# create simplified dfs to rbind 
boxplotUage <- unrepaired2019 %>% 
  as.data.frame() %>% 
  transmute(Class = recode(Class, 
                           "1" = paste0("Class 1 Leaks",
                                        " (n = ", 
                                        formatC(nrow(filter(., Class == "1")),big.mark = ","),")"),
                           "2" = paste0("Class 2 Leaks",
                                        " (n = ", 
                                        formatC(nrow(filter(., Class == "2")),big.mark = ","),")"),
                           "3" = paste0("Class 3 Leaks",
                                        " (n = ", 
                                        formatC(nrow(filter(., Class == "3")),big.mark = ","),")")),
            Utility = recode(Utility, 
                          "National Grid - Boston Gas" = "National Grid",
                          "National Grid - Colonial Gas" = "National Grid",
                          "Fitchburg" = "Unitil/Fitchburg Gas"),
            LeakAgeDays = LeakAgeDays)

# create "all leak class", assemble df and generate faceted graph
boxplotUage %>% 
  mutate(Class = paste0("All Leaks", " (n = ", 
                                        formatC(nrow(.),
                                                big.mark = ","),")")) %>% 
  rbind(., boxplotUage) %>% 
  ggplot(aes(x = reorder(Utility, LeakAgeDays, FUN = median, 
                         na.rm = TRUE), 
             y = LeakAgeDays, fill = Utility)) +
  geom_boxplot(outlier.size = 0.9) + 
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position="none") +
  ggtitle("Ages (days) of Unrepaired Leaks by Utility in 2019") +
  xlab("") + ylab("Days since leak reported to Dec 31, 2019") + 
  scale_y_continuous(labels = function(x) format(x, big.mark = ",")) +
  facet_wrap("Class")

ggsave(filename = "Images/pub/boxplotLeaksByUtilityAge.jpeg", dpi = 1000)
```


```{r figDotAgeU, fig.align='center', fig.cap="Relative exposures to population-weighted mean age of unrepaired gas leaks in 2019 by utility at Census Block Group, Census Tract, and municipality scales"}
# Set up for dot graph to show all three scales by utility
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_bg %>% 
  mutate(Group = factor(Group, levels = group_orderBG)) %>% 
  pivot_longer(c(wLeakAgeDaysAvgRRBG, wLeakAgeDaysAvgRRCG, wLeakAgeDaysAvgRREV, 
                 wLeakAgeDaysAvgRRFG, wLeakAgeDaysAvgRRLU, wLeakAgeDaysAvgRRNG), 
               names_to = "Utility", values_to = "leakDensityBG") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races",
                       "MA_MINORITY17", "MA_INCOME17")) %>% 
  drop_na(leakDensityBG) %>%
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income"),
         Utility = recode(Utility, "wLeakAgeDaysAvgRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityBG, Utility),
         Status = case_when(
           leakDensityBG == 1 ~ "Pop Mean",
           leakDensityBG > 1 ~ "Greater Exposure",
           leakDensityBG < 1 ~ "Lower Exposure"
         ),
         Scale = "Block Group") %>% 
  select(Group, leakDensityBG, Utility, Status, Scale)

# Extract rows with EJ groups to rbind with other dfs
EJrowsTR <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Tract") %>% 
  rename(leakDensityTR = leakDensityBG)

EJrowsCO <- ppLeakDensityJoinedU_bg2 %>% 
  filter(str_detect(Group, "MA ")) %>% 
  mutate(leakDensityBG = NA,
         Scale = "Municipality") %>% 
  rename(leakDensityCO = leakDensityBG)

# repeat for tracts
ppLeakDensityJoinedU_tr2 <- ppLeakDensityJoinedU_tr %>% 
  mutate(Group = factor(Group, levels = group_orderTR)) %>%
  pivot_longer(c(wLeakAgeDaysAvgRRBG, wLeakAgeDaysAvgRRCG, wLeakAgeDaysAvgRREV, 
                 wLeakAgeDaysAvgRRFG, wLeakAgeDaysAvgRRLU, wLeakAgeDaysAvgRRNG), 
               names_to = "Utility", values_to = "leakDensityTR") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityTR) %>%
  mutate(Utility = recode(Utility, "wLeakAgeDaysAvgRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityTR, Utility),
         Status = case_when(
           leakDensityTR == 1 ~ "Pop Mean",
           leakDensityTR > 1 ~ "Greater Exposure",
           leakDensityTR < 1 ~ "Lower Exposure"
         ),
         Scale = "Tract") %>% 
  select(Group, leakDensityTR, Utility, Status, Scale) %>% 
  rbind(EJrowsTR)

# repeat for cosubs
ppLeakDensityJoinedU_co2 <- ppLeakDensityJoinedU_co %>% 
  mutate(Group = factor(Group, levels = group_orderCO)) %>%
  pivot_longer(c(wLeakAgeDaysAvgRRBG, wLeakAgeDaysAvgRRCG, wLeakAgeDaysAvgRREV, 
                 wLeakAgeDaysAvgRRFG, wLeakAgeDaysAvgRRLU, wLeakAgeDaysAvgRRNG), 
               names_to = "Utility", values_to = "leakDensityCO") %>% 
  filter(!Group %in% c("Native American", "Other race", 
                       "Native Pacific Islander", "Two or more races")) %>% 
  drop_na(leakDensityCO) %>%
  mutate(Utility = recode(Utility, "wLeakAgeDaysAvgRRBG" = paste0("Berkshire Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Berkshire Gas")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRCG" = paste0("Columbia Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Columbia Gas")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRREV" = paste0("Eversource Energy"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Eversource")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRLU" = paste0("Liberty Utilities"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Liberty Utilities")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRNG" = paste0("National Grid"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "National Grid")),big.mark = ","),")"),
                          "wLeakAgeDaysAvgRRFG" = paste0("Unitil/Fitchburg Gas"," (n = ",formatC(nrow(filter(unrepaired2019,Utility == "Fitchburg")),big.mark = ","),")")),
         # Group = reorder_within(Group, leakDensityCO, Utility),
         Status = case_when(
           leakDensityCO == 1 ~ "Pop Mean",
           leakDensityCO > 1 ~ "Greater Exposure",
           leakDensityCO < 1 ~ "Lower Exposure"
         ),
         Scale = "Municipality") %>% 
  select(Group, leakDensityCO, Utility, Status, Scale) %>% 
  rbind(EJrowsCO)

# Extract rows with Disabled and HBurdened to rbind with BG df
ppLeakDensityJoinedU_bg2 <- ppLeakDensityJoinedU_co2 %>% 
  filter(str_detect(Group, "Disabled|Housing")) %>% 
  mutate(leakDensityCO = NA,
         Scale = "Block Group") %>% 
  rename(leakDensityBG = leakDensityCO) %>% 
  rbind(ppLeakDensityJoinedU_bg2,.)

# join them all together
ppLeakDensityJoinedU_all <- list(ppLeakDensityJoinedU_tr2,
                                 ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) mutate(x, UtilityGroup = paste(Utility,Group))) %>%
  reduce(., full_join, by = "UtilityGroup") %>%  
  rowwise() %>% 
  mutate(min = min(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         med = median(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         max = max(leakDensityBG, leakDensityTR, leakDensityCO, na.rm = T),
         Status = case_when(
           min > 1 & max > 1 ~ "Greater Exposure",
           min < 1 & max < 1 ~ "Lower Exposure",
           min < 1 & max > 1 ~ "Pop Mean",
           min > 1 & max < 1 ~ "Pop Mean")) %>% 
  select(-(ends_with(".x")|ends_with(".y")))

# create an alternate df_all with rbind so that we can map aesthetic to shape
ppLeakDensityJoinedU_all2 <- list(ppLeakDensityJoinedU_tr2,
                                  ppLeakDensityJoinedU_bg2,
                             ppLeakDensityJoinedU_co2) %>% 
  lapply(., function(x) rename(x,leakDensity = starts_with("leakDensity"))) %>% 
  reduce(., rbind)

# create the figure with shapes for each scale
cols <- c("#F7A35C", "#7CB5EC", "gray30")
shps <- c(0,2,1)

ppLeakDensityJoinedU_all2 %>% 
  mutate(Scale = factor(Scale, levels = c("Block Group", "Tract", 
                                          "Municipality"))) %>% 
  ggplot(aes(x = Group, y = leakDensity, color = Status, shape = Scale)) +
  geom_point() +
  geom_hline(yintercept = 1, color = "gray30") +
  scale_color_manual(values = cols, guide = "legend") +
  scale_shape_manual(values = shps, guide = "legend") +
  geom_linerange(data = ppLeakDensityJoinedU_all,
                 aes(x = Group, y = med, ymin = min, ymax = max), size = 0.3) +
  coord_flip() + 
  # scale_x_reordered() +
  theme_minimal(base_size = 6) +
  theme(legend.title=element_blank(), legend.position = "top",
        plot.caption = element_text(hjust = 0)) +
  facet_wrap(~ Utility, drop = FALSE) +
  # facet_wrap(~ Utility, scales = "free_y", drop = FALSE) +
  labs(x = NULL, 
       y = "Ratio of group population-weighted mean leak age (days) to general population-weighted mean",
       title = "Relative Exposure to Mean Age of Unrepaired Leaks in 2019 across Massachusetts by Utility", caption = "Shapes indicate scale or unit of analysis. Colors indicate if exposure is above (orange), below (blue), or equal to (gray) general population. Colors of horizontal bars through shapes\nindicate if exposure is consistently above (orange) or below (blue), or if it is mixed (gray). Environmental Justice communities ('MA') only evaluated at Block Group scale.\nDisabled Adults and Housing Burdened only evaluated at Tract and Municipality scales.")

ggsave("Images/pub/DotAgeUtilityRE.jpeg", dpi = 1000)
```


# Discussion

This analysis of natural gas leaks across Massachusetts shows that there are inequities in the geographic distribution of these leaks, as well as in how quickly they have been repaired. In general, People of Color (especially people who identify racially as Asian or Black), limited English speaking households, lower income persons, renters, and adults with lower levels of education live in neighborhoods or areas with higher leak densities, even when controlling for housing density. Moreover, these same people and places experience slower repair times, and unrepaired leaks that are significantly older than average. By contrast, people who identify racially as White and adults over 64 more often live in places with below average leak densities, and where leak repairs are completed more quickly – even within the same utility region. Whether or not these leaks represent inequitable health and safety hazards, which this analysis does not address, these inequities are significant because they reveal a procedural inequity in how leaks are being addressed. Different kinds of communities are being treated differently primarily along the axes of race, English-speaking ability, and class. This differential experience and treatment constitute an environmental injustice regardless of whether these leaks represent a health or safety hazard. 

The inequities revealed by this analysis suggest different mechanisms by which inequity occurs or is perpetrated. The metropolitan regions of Massachusetts are amongst the most racially segregated in the nation [@UMich_2010]. This segregation is a consequence of overtly discriminatory historic policies and practices, such as redlining [@Rothstein_2017], as well as ongoing and active housing discrimination in Massachusetts [@Campen_2018; @Langowski_Berman_2020], and maintained by the inertia of generational wealth differences and other structural inequalities that are the legacy of these policies and behaviors [@Killewald_Bryan_2018; @Park_Quercia_2020]. Housing discrimination and wealth inequalities have created a segregated residential geography by pushing People of Color, immigrants, and lower income groups into denser, older urban cores, as well as into cheaper and less desirable neighborhoods and housing. It has also kept them out of wealthier, more suburban, and whiter communities [@Modestino_Ziegler_2019]. The older, urban communities are also where the oldest and densest infrastructure is located, including natural gas service lines, which is where we should expect to find more gas leaks. 

But residential geography and infrastructure age alone do not explain the gas leaks inequities. The analysis also shows that unrepaired gas leaks are left unrepaired longer for the same groups that also experience higher leak densities. By contrast, leaks for Whites and those over 64 are repaired more quickly. The difference in leak age occurs primarily amongst Class 3 (i.e. less hazardous) leaks, which account for the bulk of unrepaired leaks. Although recent research has revealed that a small subset (~7%) of these leaks likely account for 50% of total distribution gas leak volume, utilities have significant discretion over these “non-hazardous” leaks. Systematic differences in the duration over which gas leaks remain unrepaired cannot simply be attributed to geographic happenstance, but may instead be related to the decision-making processes of the utilities themselves, as well as the policies (or lack thereof) that enable this unequal treatment. Beyond regulatory requirements by leak class or grade, it is less clear how utilities prioritize repairs.

This analysis stands in contrast to previous work on socioeconomic differences and gas utility service. In their analysis of the performance of natural gas utilities across the USA, @Scott_Scott_2019 did not find a clear relationship between gas service problems and the socio-economic characteristics of the utility territories. They speculated that this lack of apparent differences might be due to the scale of their analysis, which averaged socioeconomic conditions across utility territories, and thus potentially obscured smaller scale variations within those territories. The present analysis confirms those suspicions and shows that the scale of analysis does matter. Generally, differences in leak exposure between population groups are more apparent at smaller scales (i.e., smaller areas), and these differences become attenuated at larger scales. In some cases, the patterns of difference may be reordered or entirely reversed when performing the same analysis at different scales. The choice of the scale or unit of analysis should match the question and phenomenon investigated. In this case, because residential segregation occurs at sub-utility and even sub-municipal scales, it is necessary to employ units of analysis that most closely match the geography of that underlying phenomenon. 

Demographic aggregation or disaggregation affect analytic outcomes. Marginalized populations are not a homogenous group and their experiences are not interchangeable. While People of Color are generally the most impacted across the state, this impact varies by specific racial or ethnic group and by geography. Asians, Blacks, and limited English speaking households are most burdened in general. By contrast, Hispanics are not generally among the most burdened for the state as a whole, but they are in at least two utility territories. And while it is widely recognized that People of Color in the USA are disproportionately lower income, race or ethnicity and metrics of relative wealth (e.g., low income, renter vs owner) are not substitutable. Indeed, this analysis shows that race, ethnicity, and English language ability are the leading indicators of gas leaks exposure far and above indicators of wealth. Finally, in at least one utility territory, the opposite holds true: education and metrics of wealth are the leading indicators for disproportionate leak exposure, while race, ethnicity, and English language ability are not. In the Berkshire Gas territory, which serves western Massachusetts, it is less educated, lower income, and White residents who are the most burdened by gas leaks. The lesson here is that while social inequities are very real, pernicious, and nearly ubiquitous, the nature of those inequities may be geographically specific. 

One important caveat to this analysis is that utility-reported leaks are likely to be significant undercounts of actual gas leaks. Work by @PHILLIPS20131 and others has shown empirically and repeatedly that leaking natural gas in the distribution system beneath city streets in Massachusetts are much more frequent than would be apparent from data reported to the DPU [@Luna_Paulino_2018]. Separate investigations by HEET and @Edwards2021, as well as a state-authorized audit of natural gas utilities, have documented problematic inconsistencies in the way that utilities identify, monitor, and address leaks, as well as in record keeping. If the reported leaks by utilities represent a large sample, rather than the total population of leaks, then the present analysis is likely to be no less significant in revealing an inequitable problem that may be even larger than the available data show.


# Conclusion and Policy Implications

*Regulators and gas utilities should regularly evaluate their performance and structure their prioritization plans within an environmental justice framework.* Environmental inequities are real, and they are consequential. Along with a national reckoning over social inequality and systemic racism, there is growing expectation at both the state and federal levels to treat environmental justice seriously. Nationally, the Biden Administration has expressed strong support for environmental justice. It has reaffirmed the importance of the 1993 Executive Order on Environmental Justice [@ExecutiveOrder1994], created a new White House Environmental Advisory Council [@WHEJAC2021], and infused a focus on environmental justice  within the Department of Energy [@Clark_2021] and the Federal Energy Regulatory Commission [@FERC_2021; @LopezNickerson_2021]. Massachusetts has similarly reaffirmed the importance of environmental justice. Recent legislation codifies the definition of environmental justice communities and enhances outreach and notification for environmental permitting. It also requires new consideration of cumulative burdens for those communities and mandates that the DPU incorporate equity into its rules and regulations around energy infrastructure and services. We recommend that equity assessments and more inclusive planning processes should be a regular part of regulatory oversight and environmental compliance and enforcement around gas infrastructure and services. Regulators should explicitly require that utilities prioritize environmental justice communities in their planning, maintenance, and monitoring. These processes of decision making, and their outcomes, should be public, transparent, and subject to auditing or review.

*Regular and detailed public reporting are essential for assessments of equity and for accountability.* Massachusetts state policy requires utilities to submit standardized quarterly reports that allow for detailed monitoring and analysis of their performance and handling of gas leaks [@220CMR]. The analysis presented here shows very clear patterns of inequity, and these patterns would not have been detectable without those reporting and public disclosure requirements. Leak reporting should be improved by the use of GPS-derived coordinates to improve spatial accuracy and precision. We recommend that these transparency policies be mandated by the federal government to apply to utilities nationwide, and that this work be reproducible nationally. 

*Regulators and utilities need to anticipate and mitigate potentially regressive impacts of the energy transition.* Massachusetts has set a statewide goal of carbon neutrality by 2050 and the state Attorney General has directed natural gas utilities to report how they will adapt to these mandates [@Volcovici2021]. In the meantime, there is an aggressive movement across the state to push households to adopt renewable energy sources, to implement energy-saving practices and technologies, and to electrify their water and space heating appliances. The cumulative impact of these efforts will be an accelerating decline in household energy intensity and reliance on natural gas. The transition toward sustainable energy sources is happening most rapidly amongst wealthier households, especially homeowners [@Drehobl_Ross_2020; @Borenstein_2017; @Kwan_2012]. Without intervention, this energy transition may mean that the costs of maintaining the natural gas distribution system increasingly fall on a smaller pool of ratepayers in less wealthy households and communities [@Castigliego_2020]. Already regressive cost burdens and impacts could be worsened [@Zhou_Noonan_2019]. We recommend the active monitoring and assessments of current equity conditions to ensure a sustainable and just energy transition, including planning to prevent exacerbating existing inequities or creating new inequities in the future.

*Environmental justice communities should be prioritized for clean energy investment and transition efforts.* As this and other research show, these communities already bear a disproportionate share of the burdens of the current energy system. These are the communities that have the most to lose, and the most to gain, from a clean energy transition. Recent climate legislation in Massachusetts prioritizes residents of environmental justice communities for clean energy job training and placement, and minority- and women-owned businesses for access to clean energy-related startup opportunities and grants (@ClimateRoad Section 13(a)). These equity provisions are positive steps which we recommend be expanded to include development of clean energy infrastructure such as deep energy efficiency retrofits, microgrids, and geothermal district heating, and used as a national model. 

## Data Availability {-}
Supplementary materials showing intermediate validation of results, along with the datasets, R code, and the markdown document used to produce this manuscript, can be found at http://dx.doi.org/10.17632/bgx4yz67sh.2, an open-source online data repository hosted at Mendeley Data [@Luna_NicholasData].

## Funding {-}
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

# References {-}

<div id="refs"></div>



```{r tabUnrepairedBG, fig.align='center', include=FALSE}
# (APPENDIX) Appendix {-} 

# Appendix - Supplementary Figures and Tables



# Unrepaired leaks per square kilometer
ppLeakDensity_df_bg %>% 
  select(Group, wLeaksPerSqKm, wLeaksRR, wLeaksPerSqKmC1, wLeaksRRC1, wLeaksPerSqKmC2, wLeaksRRC2, wLeaksPerSqKmC3, wLeaksRRC3) %>% 
  mutate(Group = recode(Group, "MA_ENGLISH" = "MA Limited English HH",
                        "MA_MINORITY21" = "MA Minority",
                        "MA_INCOME21" = "MA Low Income")) %>%
  arrange(desc(wLeaksPerSqKm)) %>% 
  kable(longtable = T, booktabs = T,
        format.args = list(big.mark = ','), 
        caption = "Population-weighted mean leak density (leaks/sqkm) of unrepaired leaks in 2019 by Census Block Group", align = "r", digits = c(0,1,2,2,2,2,2,2,2), 
        col.names = c("Group","Per SqKm","RE","Per SqKm","RE",
                      "Per SqKm","RE","Per SqKm","RE")) %>% 
  add_header_above(., c(" ", "All Unrepaired Leaks" = 2, "Class 1 Leaks" = 2, 
                        "Class 2 Leaks" = 2, "Class 3 Leaks" = 2)) %>% 
  kable_styling(latex_options = c("repeat_header")) %>% 
  add_footnote(., "RE = Relative Exposure or ratio of group leak density to leak density of total population, total households, or total occupied housing units", notation = "none")
```