EPA Case Study

We are going apply some of the techniques we learned in Exploratory Data Analysis to study air pollution data, specifically particulate matter (we’ll call it pm25 sometimes), collected by the U.S. Environmental Protection Agency. This website https://www.health.ny.gov/environmental/indoors/air/pmq_a.htm from New York State offers some basic information on this topic if you’re interested.

Particulate matter (less than 2.5 microns in diameter) is a fancy name for dust, and breathing in dust might pose health hazards to the population. We’ll study data from two years, 1999 (when monitoring of particulate matter started) and 2012. Our goal is to see if there’s been a noticeable decline in this type of air pollution between these two years.

Reading in the data

There are 2 large files to read in using the R command read_delim. We’ll store the 1999 data in the array pm0 and let R print the tibble to see its dimensions.

library(tidyverse)

package 㤼㸱tidyverse㤼㸲 was built under R version 3.6.3Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.2.1     v purrr   0.3.3
v tibble  2.1.3     v dplyr   0.8.3
v tidyr   1.0.0     v stringr 1.4.0
v readr   1.3.1     v forcats 0.4.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

pm0 <- read_delim("RD_501_88101_1999-0.txt", comment = "#", col_names = FALSE, delim = "|", na = "")

Parsed with column specification:
cols(
  .default = col_logical(),
  X1 = col_character(),
  X2 = col_character(),
  X3 = col_character(),
  X4 = col_character(),
  X5 = col_character(),
  X6 = col_double(),
  X7 = col_double(),
  X8 = col_double(),
  X9 = col_double(),
  X10 = col_double(),
  X11 = col_double(),
  X12 = col_time(format = ""),
  X13 = col_double(),
  X14 = col_character(),
  X15 = col_double(),
  X17 = col_character()
)
See spec(...) for full column specifications.

pm0

We see that pm0 has over 117,000 lines, each containing 5 columns. In the original file, at the EPA website, each row had 28 columns, but since we’re only using a few of these, we’ve created and read in a somewhat smaller file.

We can also see there’s some missing data, but we won’t worry about that now. The column names, X1, X2, etc., are not informative. However, we know that the first line of the original file (a comment) explained what information the columns contained.

We’ll use read_lines to read in the first line of the original file and use the functions in stringr to parse out the 28 column names of the original file.

nm<-read_lines("RD_501_88101_1999-0.txt",n_max=1)
cnm<-str_split(nm,"\\|",simplify = FALSE)[[1]]
cnm<-str_replace_all(cnm," - ","")
cnm<-str_replace_all(cnm,"# ","")
cnm<-str_replace_all(cnm," ","_")
cnm<-str_replace_all(cnm,"\\(","")
cnm<-str_replace_all(cnm,"\\)","")
cnm

 [1] "RD"                                "Action_Code"                      
 [3] "State_Code"                        "County_Code"                      
 [5] "Site_ID"                           "Parameter"                        
 [7] "POC"                               "Sample_Duration"                  
 [9] "Unit"                              "Method"                           
[11] "Date"                              "Start_Time"                       
[13] "Sample_Value"                      "Null_Data_Code"                   
[15] "Sampling_Frequency"                "Monitor_Protocol_MP_ID"           
[17] "Qualifier1"                        "Qualifier2"                       
[19] "Qualifier3"                        "Qualifier4"                       
[21] "Qualifier5"                        "Qualifier6"                       
[23] "Qualifier7"                        "Qualifier8"                       
[25] "Qualifier9"                        "Qualifier10"                      
[27] "Alternate_Method_Detectable_Limit" "Uncertainty"                      

The 28 column names are all jumbled together even though they were separated by | characters. We fixed this by reassigning the column names to the output of a call to str_split (string split) with 3 arguments. The first is cnm, the pipe symbol ‘|’ is the second (use the quotation marks), and the third is the argument simplify set to FALSE.

The variable cnm now holds a list of the column headings which is nice, but we don’t need all these. We really only need State_Code, County_Code,Site_ID, Date, and Sample_Value which are X3, X4, X5, X11, and X13, respectivly. We’ll set the names of the data set to the values of our names vector, cnm, and select those variables for our the pm0 data set.

names(pm0)<-cnm
pm0 <- select(pm0,State_Code,County_Code,Site_ID,Date,Sample_Value)
pm0

Now it’s much more clear what information each column of pm0 holds. The measurements of particulate matter (pm25) are in the column named Sample_Value.

Sample_Value itself can be treated as a numeric vector (of length 117000+) with at least the first 3 values missing. Exactly what percentage of values are missing in this vector? We’ll se the R function mean with is.na() as an argument to see what percentage of values are missing (NA).

summarize(pm0,missing=mean(is.na(Sample_Value)))

So a little over 11% of the 1.1710^{5}+ are missing. We’ll keep that in mind. Now let’s start processing the 2012 data which we stored for you in the array pm1.

We’ll repeat what we did for pm0, except a little more efficiently. First read in the file and assign column names. Then we’ll print the tibble to check the dimensions.

pm1 <- read_delim("RD_501_88101_2012-0.txt", comment = "#", col_names = FALSE, delim = "|", na = "",n_max=1000000)

Parsed with column specification:
cols(
  .default = col_logical(),
  X1 = col_character(),
  X2 = col_character(),
  X3 = col_character(),
  X4 = col_character(),
  X5 = col_character(),
  X6 = col_double(),
  X7 = col_double(),
  X8 = col_double(),
  X9 = col_double(),
  X10 = col_double(),
  X11 = col_double(),
  X12 = col_time(format = ""),
  X13 = col_double(),
  X14 = col_character(),
  X15 = col_double(),
  X17 = col_character()
)
See spec(...) for full column specifications.
20268 parsing failures.
   row col           expected actual                      file
  1082 X18 1/0/T/F/TRUE/FALSE      5 'RD_501_88101_2012-0.txt'
  1091 X18 1/0/T/F/TRUE/FALSE      5 'RD_501_88101_2012-0.txt'
353841 X27 1/0/T/F/TRUE/FALSE      2 'RD_501_88101_2012-0.txt'
353842 X27 1/0/T/F/TRUE/FALSE      2 'RD_501_88101_2012-0.txt'
353843 X27 1/0/T/F/TRUE/FALSE      2 'RD_501_88101_2012-0.txt'
...... ... .................. ...... .........................
See problems(...) for more details.

names(pm1)<-cnm
pm1 <- select(pm1,State_Code,County_Code,Site_ID,Date,Sample_Value)
pm1

Wow! nrow(pm1) rows! Particulate matter was first collected in 1999 so perhaps there weren’t as many sensors collecting data then as in 2012 when the program was more mature.

Now let’s see what percentage of values are missing in measurements of particulate matter (pm25) for 2012. As before, use the R function mean with is.na() as an argument to find out.

summarize(pm1,missing=mean(is.na(Sample_Value)))

So only `round(mean(is.na(pm1$Sample_Value)*100,1)% of the particulate matter measurements are missing. That’s about half the percentage as in 1999.

Now let’s look at summaries (using the summary command) for both datasets.

summary(pm0$Sample_Value)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   0.00    7.20   11.50   13.74   17.90  157.10   13217 

summary(pm1$Sample_Value)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 -10.00    4.20    7.90    9.25   12.00  893.70   52608 

The numbers in the vectors pm0$Sample_Value and x1pm1$Sample_Value` represent measurements taken in micrograms per cubic meter.

We see that both the median and the mean of measured particulate matter have declined from 1999 to 2012. In fact, all of the measurements, except for the maximum and missing values (Max and NA’s), have decreased. Even the Min has gone down from NA to NA! We’ll address what a negative measurment might mean a little later. Note that the Max has increased from NA in 1999 to NA in 2012. This is quite high and might reflect an error in the table or malfunctions in some monitors.

Making some exploratory plots

To make this data a bit easier to work with, I’m going to put them together. I’ll create a factor variable called Year to distinguish the measurements in 1999 from 2012. Then I’ll use the bind_rows function to sandwhich these two data frames together.

pm25<-bind_rows(pm0,pm1)
pm25<-mutate(pm25,Year=as_factor(str_sub(Date,1,4)))
pm25

Let’s make a boxplot function with 2 arguments, x0 and x1.

library(plotly)

package 㤼㸱plotly㤼㸲 was built under R version 3.6.3Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: 㤼㸱plotly㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

p<-ggplot(pm25, aes(x = Year, y = Sample_Value)) + 
  geom_boxplot(notch = TRUE) +
  stat_summary(fun.y=mean,geom="point",shape=23,size=4) # add mean points
ggplotly(p)

Removed 65825 rows containing non-finite values (stat_boxplot).Removed 65825 rows containing non-finite values (stat_summary).

Huh? Did somebody step on the boxes? It’s hard to see what’s going on here. There are so many values outside the boxes and the range of particulate matter (pm25) for 2012 is so big that the boxes are flattened. It might be more informative to call boxplot on the logs (base 2).

p<-ggplot(pm25, aes(x = Year, y = log2(Sample_Value))) + 
  geom_boxplot(notch = TRUE) +
  stat_summary(fun.y=mean,geom="point",shape=23,size=4) # add mean points
ggplotly(p)

NaNs producedNaNs producedNaNs producedRemoved 109290 rows containing non-finite values (stat_boxplot).Removed 109290 rows containing non-finite values (stat_summary).

A bonus! Not only do we get a better looking boxplot we also get some warnings from R in Red. These let us know that some values in particulate matter (pm25) for both years were “unloggable”, no doubt the 0 (Min) we saw in the summary of particulate matter (pm25) in 1999 and the negative values we saw in the Min of the summary of particulate matter (pm25) for 2012.

From the boxplot, what can you say about the data? 1: The boxes are too small to interpret 2: The range of particulate matter (pm25) in 1999 is greater than the range of particulate matter (pm25) in 2012 3: The median of particulate matter (pm25) in 2012 is less than the median of particulate matter (pm25) in 1999 4: The mean of particulate matter (pm25) in 2012 is less than the mean of particulate matter (pm25) in 1999

Let’s return to the question of the negative values in particulate matter (pm25) for 2012. Let’s count how many negative values there are. We’ll do this in a few steps.

First, group the data by year. Then, form the vector negative by assigning to it the boolean Sample_Value<0. Now run the R command summarize. We can use the sum function and set the na.rm argument equal to TRUE. This tells sum to ignore the missing values in our logic.

pmyr<-group_by(pm25,Year)
summarize(pmyr,count_neg=sum(Sample_Value<0,na.rm=TRUE))

So there are over 210^{4} negative values. Sounds like a lot. Is it? Let’s run the same R command using mean to calculate a percentage missing.

summarize(pmyr,pneg=mean(Sample_Value<0,na.rm=TRUE))

We see that just 2.2% of the particulate matter (pm25) for 2012 values are negative. Perhaps that’s a small enough percentage that we can ignore them. Before we ignore them, though, let’s see if they occur during certain times of the year.

But Date is a very long set of integers, and the format of the entries is hard to read. There’s no separation between the year, month, and day. Let’s reassign the output of a call to ymd.

library(lubridate)

Attaching package: 㤼㸱lubridate㤼㸲

The following object is masked from 㤼㸱package:base㤼㸲:

    date

pm25<-mutate(pm25,Date=ymd(Date))
pm25

Now Date is in a nicer format.

Let’s plot the months when the particulate matter measurements are negative.

pm_2k12_neg<-filter(pm25,Year=='2012',Sample_Value<0)
p<-ggplot(pm_2k12_neg, aes(x=as_factor(month(Date)))) +
  geom_bar(fill="#bcbddc") +
  xlab("Month")
ggplotly(p)

We see the bulk of the negative measurements were taken in the winter months, with a spike in May. Not many of these negative measurements occurred in summer months. We can take a guess that because particulate measures tend to be low in winter and high in summer, coupled with the fact that higher densities are easier to measure, that measurement errors occurred when the values were low. For now we’ll attribute these negative measurements to errors. Also, since they account for only 2% of the 2012 data, we’ll ignore them.

New York State

Now we’ll change focus a bit and instead of looking at all the monitors throughout the country and the data they recorded, we’ll try to find one monitor that was taking measurements in both 1999 and 2012. This will allow us to control for different geographical and environmental variables that might have affected air quality in different areas. We’ll narrow our search and look just at monitors in New York State.

First we’ll subset off the New York monitor identification data for 1999 and 2012 into 2 vectors, site0 and site1. Look at the structure of site0 now with the R command str.

pm25<-mutate(pm25,Site=as_factor(str_c(County_Code,Site_ID)))
site1999<-filter(pm25,State_Code=='36',Year=='1999')
site2012<-filter(pm25,State_Code=='36',Year=='2012')
site1999

site2012

We see that the IDs of monitors in New York State in 1999 is a set of strings, each of which has the form CCCSSSS. We created these from the county codes (the C portion of the string) and the monitor IDs (the S portion). We can use the intersect command with to examine which sites have measurements in both years.

both <- dplyr::intersect(site1999$Site, site2012$Site)
both

 [1] "0010005" "0010012" "0050080" "0130011" "0290005" "0310003" "0632008" "0671015"
 [9] "0850055" "1010003"

We see that 10 monitors in New York State were active in both 1999 and 2012.

Now let’s see how many measurements each of the 10 New York monitors that were active in both 1999 and 2012 took in those years. We’ll create 2 subsets (one for each year), one of pm0 and the other of pm1.

The subsets will filter for 2 characteristics. The first is State.Code equal to 36 (the code for New York), and the second is that the county.site (the component we added) is in the vector both.

First create the variable cnt0 by assigning to it the output of the R command subset, called with 2 arguments. The first is pm0, and the second is a boolean with the 2 conditions we just mentioned. Recall that the testing for equality in a boolean requires ==, intersection of 2 boolean conditions is denoted by & and membership by %in%.

Take a look at both now.

site<-bind_rows(site1999, site2012)
site <- filter(site, Site %in% both)
site

Now run the split and map trick we saw in the Iteration lecture. This will split Site into several data frames according to county.site (that is, monitor IDs) and tell us how many measurements each monitor recorded.

site %>%
  group_by(Year,Site) %>%
  summarize(count=n(),
            avg_pm25=mean(Sample_Value,na.rm=TRUE))

From the output above, which monitor is the only one whose number of measurements increased from 1999 to 2012? 1: 0632008 2: 1010003 3: 0290005 4: 0850055 5: None

We want to examine a monitor with a reasonable number of measurements so let’s look at the monitor with ID 0632008. Create a variable pm0sub which is the subset of cnt0 (this contains just New York data) which has County.Code equal to 63 and Site.ID 2008.

site.sub <- filter(site, Site=='0632008')
site.sub

From the output of the 2 calls to sapply, how many rows will site.sub have for 1999? 1: 183 2: 29 3: 122 4: 30

Now we’d like to compare the pm25 measurements of this particular monitor (0632008) for the 2 years. First, create the vector x0sub by assigning to it the Sample.Value component of pm0sub.

g<-g<-ggplot(site.sub,aes(Year,Sample_Value)) +
  geom_violin() +
  stat_summary(fun.y=median, geom="crossbar", size=2, color="#a6bddb")
ggplotly(g)

Removed 35 rows containing non-finite values (stat_ydensity).Removed 35 rows containing non-finite values (stat_summary).

Which median is larger - the one for 1999 or the one for 2012? 1: 2012 2: 1999

We’ve created a violin plot showing the PM 2.5 values side by side with the same range of values on the y axis. The improvement in the medians between 1999 and 2012 is now clear. Also notice that in 2012 there are no big values (above 15). This shows that not only is there a chronic improvement in air quality, but also there are fewer days with severe pollution.

The last avenue of this data we’ll explore (and we’ll do it quickly) concerns a comparison of all the states’ mean pollution levels. This is important because the states are responsible for implementing the regulations set at the federal level by the EPA.

Let’s first gather the mean (average measurement) for each state in both years.

by_state<-pm25 %>%
  group_by(State_Code,Year) %>%
  summarise(avgpm25=mean(Sample_Value,na.rm=TRUE)) %>%
  pivot_wider(names_from = Year, values_from = avgpm25) %>%
  mutate(increase=`1999`<`2012`)
by_state

Each row of mrg has 4 entries - a state identified by number, a state mean for 1999, a state mean for 2012, and an indicator for whether or not the state saw an increase in PM 2.5 levels between 1999 and 2012.

Now we’ll plot the data to see how the state means changed between the 2 years. I’m going to use a parallel coordinate plot using ggally

library(GGally)

package 㤼㸱GGally㤼㸲 was built under R version 3.6.3Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

Attaching package: 㤼㸱GGally㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    nasa

p <- ggparcoord(data = by_state, columns = 2:3,groupColumn =4,showPoints = TRUE,title ="4 States had increased PM 2.5 levels between 1999 and 2012",scale="globalminmax")+
  scale_colour_manual(values = c('#bdbdbd','#de2d26')) + 
  theme(legend.position = "none") +
  ylab("PM 2.5 concentration") +
  xlab("Year")
ggplotly(p)

We see a column of points at 1999 which represent the 1999 state means. We see a shorter column of points at 2012. We see from the plot that the vast majority of states have indeed improved their particulate matter counts so the general trend is downward. There are a few exceptions. (The topmost point in the 1999 column is actually two points that had very close measurements.)

For fun, let’s see which states had higher means in 2012 than in 1999. Just use the mrg[mrg\(mean.x < mrg\)mean.y, ] notation to find the rows of mrg with this particulate property.

filter(by_state,increase)

Only 4 states had worse pollution averages, and 2 of these had means that were very close. If you want to see which states (15, 31, 35, and 40) these are, you can check out this website https://www.epa.gov/enviro/state-fips-code-listing to decode the state codes.

This concludes our comparison of air pollution data from two years in different ways. First, we looked at measures of the entire set of monitors, then we compared the two measures from a particular monitor, and finally, we looked at the mean measures of the individual states.

Congratulations! We hope you enjoyed this particulate lesson.