What Is PM2.5 and Why Does It Matter?

Fine particulate matter, or PM2.5, refers to airborne particles smaller than 2.5 micrometers in diameter — roughly 1/30th the width of a human hair. The “2.5” is the upper size limit in micrometers (µm); one micrometer is one millionth of a meter. Because these particles are so small, they stay suspended in the air for days, can travel hundreds of kilometers on the wind, and — most importantly — bypass the body’s normal defenses. When inhaled, they penetrate deep into the lungs and can cross into the bloodstream.

Long-term exposure is linked to respiratory disease, cardiovascular disease, and premature death. Short-term spikes — a smoke-filled day from a nearby wildfire, for example — can trigger asthma attacks and hospital admissions within hours. PM2.5 is not uniformly distributed: concentrations are highest near major roads, industrial facilities, and agricultural burning, which means that communities living in those places bear a disproportionate health burden.

Where Does PM2.5 Come From?

Sources fall into two broad categories. Primary PM2.5 is emitted directly as particles: diesel exhaust, wood smoke, dust from unpaved roads, and agricultural operations. Secondary PM2.5 forms in the atmosphere when gases — sulfur dioxide from power plants, nitrogen oxides from vehicles — react chemically and condense into particles.

Concentrations are also shaped by meteorology. Temperature inversions trap pollution near the surface on still winter mornings; onshore sea breezes and afternoon convection can flush coastal basins; precipitation scavenges particles from the air. These processes produce strong seasonal cycles that vary by region.

What This Analysis Does

Because PM2.5 patterns depend so heavily on local geography and emission sources, each student in this class analyzes a different county so we can compare results across regions. Using thirty years of EPA monitoring data, you will:

• Describe the monitoring network in your county and the stations that have long enough records to support analysis.
• Visualize the full daily record to identify long-term trends, seasonal patterns, and episodic spikes.
• Quantify monthly averages and their uncertainty using 95% confidence intervals.
• Evaluate how often concentrations exceeded the federal health standards, and whether air quality has improved over the past three decades.

Regulatory Standards

The EPA sets two National Ambient Air Quality Standards (NAAQS) for PM2.5. Understanding both is essential for interpreting your results.

The 24-hour standard is 35 ug/m3. No more than one exceedance per year is permitted on a three-year rolling average. It protects against acute exposure — a single smoky day from a nearby wildfire.

The annual standard is 12 ug/m3. The three-year mean of annual averages must stay below this value. It protects against chronic exposure: living near a freeway for years.

All time-series plots in this document carry a reference line at 35 ug/m3 because our data are daily values.

What Your Submission Should Include

For your own document, include:

• Thirty years of daily PM2.5 data downloaded and loaded for your chosen county, with a brief description of how many stations were retained and why any were removed.
• A summary statistics table (median, mean, SD, 95% CI by month) for at least one station, with a sentence interpreting the seasonal pattern.
• A scatter plot of all daily values with the 24-hour NAAQS annotated, and a monthly mean +/- 95% CI plot for your focus station.
• A discussion of 2–4 paragraphs addressing regional meteorology, emissions sources, exceedances, and environmental justice.

Project Directory Setup

You are working on a shared RStudio Server. Each student has their own home directory (/home/yourusername/). Keep your work there — do not write files to shared system directories.

Create your project folders using the Files pane (lower-right panel in RStudio):

1. In the Files pane, make sure you are in your home directory — click the house icon (🏠) if you are not sure.
2. Click New Folder, type EJnPi, and click OK.
3. Double-click into EJnPi to open it.
4. Click New Folder again, type data, and click OK.

Then:

1. Create a new .Rmd file: File → New File → R Markdown, give it any title, and click OK.
2. Immediately save it: File → Save As, navigate into the EJnPi folder, and name it Analyzing_EPA_PM25_YourLastName.Rmd.
3. At the top of your .Rmd, add this line so knitr always runs in the right place:

knitr::opts_knit$set(root.dir = "~/EJnPi")

Your folder structure before the first knit:

~/EJnPi/
  Analyzing_EPA_PM25_YourLastName.Rmd
  data/              <- empty; the download loop will populate it

Course Functions

This assignment uses R functions written for this course that handle downloading and cleaning the EPA data, so you can focus on the analysis. They are maintained in the course GitHub repository.

Browse the code: https://github.com/marclos/EJnPi/tree/master/docs

You do not need to download or copy anything. Your Rmd loads them automatically with this line (already in the setup section below):

source("https://raw.githubusercontent.com/marclos/EJnPi/master/docs/EJnPi_functions.R")

If you are curious how the functions work, click the link above, open EJnPi_functions.R, and read the comments. You are not expected to modify them.

Install Required Packages

On the server, packages install to your personal library (~/R/). Run this once in the RStudio Console — not inside your .Rmd — before knitting for the first time:

install.packages(c("dplyr", "lubridate", "xtable"))

Find Your FIPS Codes

EPA data files are indexed by state and county FIPS codes — standardized numeric identifiers used across all federal datasets. You need a 2-digit state code and a 3-digit county code, both zero-padded (California = "06", not "6").

1. Go to: https://www.census.gov/library/reference/code-lists/ansi.html#cou
2. Click your state to open a CSV lookup table.
3. Find your county in COUNTYNAME; note the value in COUNTYFP.
4. The state code (STATEFP) is in the first column.

Load Libraries

Always load packages at the very top of your script so dependencies are explicit and any missing-package errors appear before your analysis starts. Also set the working directory in your setup chunk so knitr resolves file paths correctly on the server:

knitr::opts_knit$set(root.dir = "~/EJnPi")   # adjust to your actual folder name

library(dplyr)       # group_by(), summarise(), mutate(), filter()
library(lubridate)   # ymd(), month() -- date parsing and extraction

Load EPA Data

The EPA Air Quality System (AQS) is the federal database of air pollution measurements from thousands of monitoring stations across the United States. Stations report daily average concentrations; the data are freely available as annual CSV files at a stable public URL — no account or API key required.

The variable we analyze is PM2.5 FRM/FEM mass (EPA pollutant code 88101), measured in micrograms per cubic meter (µg/m³). Each annual file covers the entire country; the helper function filters it to your county using the FIPS codes you configured above and caches the result locally so you only download once.

A helper function handles the download and caching. It is defined below in a hidden chunk — you do not need to read or modify it. On the first knit it downloads thirty years of data (allow 20–30 minutes on first knit). On every subsequent knit it finds the cached file and loads it in seconds.

The helper functions live in a shared file on GitHub. The chunk below sources them directly — no download or installation needed. If the server has no outbound internet access, ask your instructor for the local path.

source("https://raw.githubusercontent.com/marclos/EJnPi/master/docs/EJnPi_functions.R")

get_epa_pm25() handles downloading and caching; internally it calls prepare_epa_pm25() to select and rename columns. Call just one line:

end_year   <- as.integer(format(Sys.Date(), "%Y")) - 1
start_year <- end_year - 29   # 30 years

# I name my dataframe SB for San Bernardino.
# In your document, replace SB with your chosen abbreviation everywhere.
SB <- get_epa_pm25(
  state_fips  = state_fips,
  county_fips = county_fips,
  start_year  = start_year,
  end_year    = end_year
)

## Loading cached data: data/06071_PM25_daily.rds

Inspect the Data

Always look at a new dataframe before doing anything with it. The function returns five columns. Confirm that Date is class Date (not character) and that PM2.5.AQI.Value is numeric. If either looks wrong, plots and statistics will silently produce incorrect results.

str(SB)

## 'data.frame':    21998 obs. of  5 variables:
##  $ Site.ID        : chr  "06-071-0014" "06-071-0014" "06-071-0014" "06-071-0014" ...
##  $ Site.Name      : chr  "Victorville-Amargosa Road" "Victorville-Amargosa Road" "Victorville-Amargosa Road" "Victorville-Amargosa Road" ...
##  $ Date           : Date, format: "1999-01-03" "1999-01-09" ...
##  $ PM2.5.AQI.Value: int  62 60 58 57 52 65 56 55 66 34 ...
##  $ PM2.5.Conc     : num  15.2 13.9 12.7 12.1 9.9 16.7 11.7 11.1 17.4 6.2 ...

A quick look at the first few rows confirms the column names and values are sensible:

head(SB)

##       Site.ID                 Site.Name       Date PM2.5.AQI.Value PM2.5.Conc
## 1 06-071-0014 Victorville-Amargosa Road 1999-01-03              62       15.2
## 2 06-071-0014 Victorville-Amargosa Road 1999-01-09              60       13.9
## 3 06-071-0014 Victorville-Amargosa Road 1999-01-15              58       12.7
## 4 06-071-0014 Victorville-Amargosa Road 1999-01-21              57       12.1
## 5 06-071-0014 Victorville-Amargosa Road 1999-01-27              52        9.9
## 6 06-071-0014 Victorville-Amargosa Road 1999-01-30              65       16.7

Assign Decade Blocks

With 30 years of data we can ask whether air quality has changed over time. We split the record into three 10-year blocks and store the label as a factor so plots order the decades correctly.

# Extract the 4-digit year, round down to the decade start, build a label
SB$Year   <- as.integer(format(SB$Date, "%Y"))
SB$Decade <- factor(
  paste0(floor(SB$Year / 10) * 10, "s"),   # e.g. 2010 -> "2010s"
  levels = sort(unique(paste0(floor(SB$Year / 10) * 10, "s")))
)

table(SB$Decade)   # confirm three roughly equal groups

## 
## 1990s 2000s 2010s 2020s   NAs 
##   532  1496  9345  3371  7254

Remove QA Collocated Monitors

Some monitoring sites are quality-assurance collocations — duplicate instruments installed next to a primary monitor to audit its performance. Including them double-counts measurements from those locations. They are identified by site numbers in the 9000-range (fourth segment starting with 9).

SB <- subset(SB, !(Site.ID %in% c("06-071-9010", "06-071-9002")))

Remove Low-Sampling Stations

A station with very few observations over thirty years cannot support reliable monthly averages. We count observations per station, print the table so you can see the full picture, then filter.

# aggregate() applies a function (length = count) to one column grouped by another
SB.N           <- aggregate(PM2.5.AQI.Value ~ Site.ID, data = SB, FUN = length)
names(SB.N)[2] <- "N"
print(SB.N)

##       Site.ID    N
## 1 06-071-0014  114
## 2 06-071-0025  914
## 3 06-071-0027 7003
## 4 06-071-0306 7410
## 5 06-071-2002 2170
## 6 06-071-8001 2286
## 7 06-071-9004 2101

# Keep only stations that contributed at least 100 daily observations
keep_ids <- SB.N$Site.ID[SB.N$N >= 100]
SB       <- subset(SB, Site.ID %in% keep_ids)

Now that the data are clean, we move to exploration. The goal of the next section is not yet to draw conclusions — it is to understand the shape of the data before computing any statistics.

Exploratory Graphics

Never skip exploratory visualization. Plotting the raw data before computing statistics helps you spot outliers, record gaps, and unusual station behavior that summary numbers alone would hide.

Station Name Lookup Table

Before plotting we create a small table mapping each Site.ID to its human-readable Site.Name. A station ID alone is not informative in a figure legend.

# unique() on two columns returns one row per distinct combination
SB.unique <- unique(SB[c("Site.ID", "Site.Name")])
SB.unique

##          Site.ID                  Site.Name
## 1    06-071-0014  Victorville-Amargosa Road
## 115  06-071-0025       Ontario Fire Station
## 211  06-071-2002                    Fontana
## 332  06-071-8001                   Big Bear
## 429  06-071-9004             San Bernardino
## 635  06-071-0306    Victorville-Park Avenue
## 4022 06-071-0027 Ontario-Route 60 Near Road

Box Plot by Station

A box plot shows the full distribution — median, interquartile range, and outliers — for every station simultaneously. The horizontal reference line at 35 ug/m3 makes it easy to see which stations routinely approach or exceed the 24-hour standard.

Here is the output your code should produce:

Distribution of daily PM2.5 AQI by monitoring station. Box = IQR; line = median; points = outliers. Red dashed line = 24-hr NAAQS (35 ug/m3).

What this chart shows. Each box represents one monitoring station’s full thirty-year record of daily PM2.5 values. The thick horizontal line inside each box is the median — half of all daily readings fell below this value, half above. The box itself spans the middle 50% of readings (the interquartile range, or IQR). Points above the whiskers are outliers: unusually high days, often driven by wildfires or pollution events.

What to look for in your data. Are some stations consistently higher than others? Do any stations have their median near or above the red dashed line? Wide boxes indicate high day-to-day variability; tall outlier stacks suggest occasional severe episodes. These differences between stations often reflect proximity to roads, industrial zones, or the geography of valley air basins.

Scaffold to adapt for your own document.

# par() sets margins and label orientation before the plot is drawn.
# Skipping it often causes station labels to be clipped — mentors see this a lot.
par(las = ___, mar = c(___, ___, ___, ___))

boxplot(PM2.5.AQI.Value ~ Site.ID, data = ___,   # your dataframe abbreviation, e.g. SB
        main = paste("PM2.5 by Station --", county),
        ylab = ___,
        xlab = ___,
        ylim = c(0, ___),
        col  = ___)

abline(h = ___, col = ___, lty = ___, lwd = ___)

legend(___, legend = ___, col = ___, lty = ___, bty = "n")

Rank Stations by Median PM2.5

Before choosing a focus station for the monthly analysis, rank all stations by their median PM2.5. This identifies both the cleanest and most-polluted monitoring locations in your county.

SB.median <- aggregate(PM2.5.AQI.Value ~ Site.ID, data = SB, FUN = median)
SB.median <- merge(SB.median, SB.unique, by = "Site.ID")
SB.median <- SB.median[order(SB.median$PM2.5.AQI.Value), ]
print(SB.median)

##       Site.ID PM2.5.AQI.Value                  Site.Name
## 6 06-071-8001              31                   Big Bear
## 4 06-071-0306              46    Victorville-Park Avenue
## 1 06-071-0014              55  Victorville-Amargosa Road
## 7 06-071-9004              55             San Bernardino
## 5 06-071-2002              56                    Fontana
## 2 06-071-0025              58       Ontario Fire Station
## 3 06-071-0027              59 Ontario-Route 60 Near Road

Time-Series: All Stations

Plotting every daily value over thirty years reveals long-term trends, year-to-year variability, and episodic spikes (often wildfire smoke). The fourth argument of rgb() sets transparency (0 = invisible, 1 = solid); using 0.25 prevents overplotting when thousands of points stack on top of each other.

Here is the output your code should produce:

Daily PM2.5 AQI for all retained stations over 30 years. Semi-transparent points reveal density. Red dashed line = 24-hr NAAQS (35 ug/m3).

What this chart shows. Every daily measurement from every retained station, plotted over thirty years. Semi-transparent points prevent the high-density periods from becoming a solid black blob — where points are darker, many readings stack on top of each other. The red dashed line marks the 24-hour NAAQS at 35 µg/m³.

What to look for. A downward drift in the cloud of points over time indicates improving air quality — a real trend seen in many US counties since the Clean Air Act was strengthened in the 1990s and 2000s. Vertical spikes above the reference line are exceedance events; if they cluster in late summer, wildfire smoke is a likely cause. Gaps in the record reveal when stations were offline.

Scaffold:

par(las = 1)

plot(___ ~ ___, data = ___,   # y ~ x; replace ___ with your dataframe abbreviation
     xlab = ___,
     ylab = ___,
     main = paste("PM2.5 --", county, "(All Stations)"),
     pch  = ___,
     cex  = ___,
     col  = rgb(0, 0, 0, ___))

abline(h = ___, col = ___, lty = ___)

Time-Series by Station

Separating stations into panels reveals whether all sites share the same seasonal pattern or whether some behave differently — a clue to localized sources or topographic effects. par(mfrow = c(nrows, ncols)) arranges plots in a grid; always reset it with par(mfrow = c(1, 1)) when done.

This section uses a for loop — one of the most powerful ideas in programming. Instead of writing seven nearly-identical plot() calls by hand, the loop runs the same block of code once for each station, swapping in a different site each time. If you had 50 stations, the code would be exactly the same length. That is the point: loops let the computer do the repetitive work so you don’t have to.

That said, I strongly recommend writing and running the code for a single station first before wrapping it in a loop. Pick one Site.ID, subset the data manually, and get a plot you are happy with. Once that works, it is straightforward to drop the same logic inside the loop. Debugging a broken loop is much harder than debugging a single plot.

Here is the output your code should produce:

Daily PM2.5 per station. Each panel is one site. Compare the magnitude and timing of peaks across locations. I used par(mfrow=c(3,3)) for San Bernardino — adjust the grid to match your number of stations.

What this chart shows. The same time-series as above, but split into one panel per station. This makes it possible to see whether all stations rise and fall together (county-wide events like regional smoke) or whether some stations behave differently from the rest (local sources, topographic effects, or data gaps at specific monitors).

What to look for. Do peaks in one station align with peaks in others? Is one station consistently higher across all decades? Does any station have a long gap where no data were collected? These observations should inform which station you choose as your focus site in the next section.

With the data explored, we now move to quantitative summaries. The goal shifts from seeing the shape of the data to measuring it precisely, with uncertainty estimates that let us compare months and decades rigorously.

Scaffold — note that a for loop iterates over a vector; inside the loop, i takes each value in turn:

site_ids <- unique(___$Site.ID)   # replace ___ with your dataframe abbreviation
n_sites  <- length(site_ids)
ncols    <- min(3, n_sites)
nrows    <- ceiling(n_sites / ncols)
colors   <- rainbow(n_sites)

par(las = 1, mfrow = c(nrows, ncols), mar = c(3, 4, 3, 1))

for (i in seq_along(site_ids)) {
  sub        <- subset(___, Site.ID == ___)
  # ^ second ___: the loop counter i moves through site_ids one at a time.
  #   On the first iteration i=1, on the second i=2, etc.
  #   So to get the current station's ID, index site_ids with i: site_ids[i]

  site_label <- ___.unique$Site.Name[___.unique$Site.ID == ___]
  # ^ last ___: same answer — site_ids[i]. You are asking:
  #   "find the row in ___.unique where Site.ID matches this iteration's station,
  #   and return its Site.Name." The answer to both blanks in this loop is the same thing.

  plot(___ ~ ___, data = sub,
       main = ___,
       xlab = "",
       ylab = "PM2.5 AQI (ug/m3)",
       pch  = 20,
       cex  = 0.5,
       ylim = c(0, ___),
       col  = colors[i])

  abline(h = ___, col = ___, lty = ___)
}

par(mfrow = c(1, 1))

Monthly Summary Statistics

Aggregating by month across all years reveals the seasonal cycle. dplyr’s group_by() + summarise() pipeline computes grouped statistics in one readable block. .groups = "drop" removes the grouping structure after summarising so the result is a plain dataframe.

# I use SB here — replace with your abbreviation in your own document
Monthly <- SB %>%
  mutate(Month = month(Date)) %>%         # integer month 1-12
  group_by(Decade, Month, Site.ID) %>%    # grouped by decade so we can compare blocks
  summarise(
    median  = median(PM2.5.AQI.Value, na.rm = TRUE),
    mean    = mean(PM2.5.AQI.Value,   na.rm = TRUE),
    sd      = sd(PM2.5.AQI.Value,     na.rm = TRUE),
    N       = n(),
    .groups = "drop"
  )

head(Monthly)

## # A tibble: 6 × 7
##   Decade Month Site.ID     median  mean    sd     N
##   <fct>  <dbl> <chr>        <dbl> <dbl> <dbl> <int>
## 1 1990s      1 06-071-0014    59   59    4.47     6
## 2 1990s      1 06-071-0025   102. 110.  33.7     10
## 3 1990s      1 06-071-2002    84   91.7 37.1     17
## 4 1990s      1 06-071-9004    67   69   28.7      5
## 5 1990s      2 06-071-0014    57   58.8 12.8      9
## 6 1990s      2 06-071-0025    81   87.3 26.3      9

Where Does 1.96 Come From?

The multiplier 1.96 is the critical value of the standard normal distribution that cuts off 2.5% in each tail, leaving 95% in the middle. The plot below shows this geometrically: the shaded regions are the rejection regions for a two-tailed test at α = 0.05.

Standard normal density curve. Red shading = rejection regions at α = 0.05 (2.5% per tail). The vertical dashed lines fall at ±1.96, the critical values that define the boundary. In R: qnorm(0.975) = 1.96.

Confidence Interval Formula

The sample mean \(\bar{x}\) estimates the true long-run monthly average, but carries uncertainty. The 95% confidence interval quantifies it:

\[CI_{95} = \bar{x} \pm 1.96 \times \frac{s}{\sqrt{N}}\]

If we repeated this sampling process many times, 95% of the resulting intervals would contain the true mean. Narrow CI = many observations have pinned the mean down precisely.

# qnorm(0.975) = 1.96 -- the z critical value shown in the plot above
Monthly$UCL95 <- Monthly$mean + qnorm(0.975) * Monthly$sd / sqrt(Monthly$N)
Monthly$LCL95 <- Monthly$mean - qnorm(0.975) * Monthly$sd / sqrt(Monthly$N)

Choose a Focus Station

For the table and CI plot we work with one station. The median-ranked station is a representative middle case — not the best or worst.

# ceiling() ensures we always land on a valid row index
focus_site <- SB.median$Site.ID[ceiling(nrow(SB.median) / 2)]
focus_site

## [1] "06-071-9004"

Summary Table

xtable converts a dataframe to formatted HTML or LaTeX. The digits vector controls decimal places column by column: position 1 is the (unused) row name; positions 2 onward correspond to each data column left to right.

Here is the table your code should produce:

Scaffold:

library(xtable)

tbl       <- Monthly[Monthly$Site.ID == focus_site, ]
tbl$Month <- month.abb[tbl$Month]   # converts 1-12 to Jan-Dec
tbl       <- tbl[, c("Month", "N", "median", "mean", "sd", "LCL95", "UCL95")]

print(
  xtable(tbl,
         caption = ___,                              # quoted caption text
         digits  = c(0, 0, ___, ___, ___, ___, ___)),# one value per column; 0 = no decimals
  include.rownames  = FALSE,
  type              = "html",
  caption.placement = "top"
)

Monthly Mean +/- 95% CI by Decade

With 30 years of data we can ask: has PM2.5 changed over time? We plot three lines — one per decade — on the same axes. Lines that shift downward over time indicate improving air quality; lines that converge or diverge reveal whether seasonal patterns have changed. Error bars are 95% CIs on the monthly mean, so non-overlapping bars between decades indicate a statistically meaningful difference.

Here is the output your code should produce:

Monthly mean PM2.5 AQI +/- 95% CI for the focus station, split by decade. Each line is one 10-year block. Non-overlapping error bars between decades indicate a statistically meaningful shift in that month’s typical concentration. Red dashed line = 24-hr NAAQS (35 ug/m3).

What this chart shows. Each line represents one decade’s monthly average PM2.5 concentration at the focus station, with error bars showing the 95% confidence interval around each mean. Darker lines are older decades; lighter lines are more recent.

What to look for. Lines shifting downward from dark to light indicate improving air quality over time — a statistically meaningful improvement if the error bars do not overlap between decades. The shape of each line reveals the seasonal cycle: which months are typically worst and which are cleanest, and whether that pattern has changed across decades. A seasonal peak in late summer often signals wildfire smoke; a winter peak often reflects temperature inversions trapping vehicle emissions near the surface.

Scaffold — the key new ideas here are plot(NULL, ...) to set up empty axes before drawing anything, and the offset jitter so overlapping error bars stay readable.

focus   <- Monthly[Monthly$Site.ID == focus_site, ]
decades <- levels(focus$Decade)
n_dec   <- length(decades)

pal      <- colorRampPalette(c(___, ___, ___))(n_dec)
ylim_top <- max(focus$UCL95, 35, na.rm = TRUE) * 1.15

par(las = ___, mar = c(___, ___, ___, ___))
plot(NULL,
     xlim = c(0.5, 12.5), ylim = c(0, ylim_top),
     xaxt = "n",
     xlab = ___, ylab = ___, main = ___)

axis(1, at = 1:12, labels = month.abb)
abline(h = ___, col = ___, lty = ___, lwd = ___)

offsets <- seq(-0.15, 0.15, length.out = n_dec)

for (i in seq_along(decades)) {
  d  <- focus[focus$Decade == decades[i], ]
  d  <- d[order(d$Month), ]
  mx <- d$Month + offsets[i]

  lines(mx, d$___, col = pal[i], lwd = ___)
  points(mx, d$___, col = pal[i], pch = ___)

  arrows(mx, d$LCL95, mx, d$UCL95,
         angle = 90, code = 3, length = ___, col = pal[i], lwd = ___)
}

legend(___, legend = c(decades, ___),
       col = c(pal, ___), lty = ___, pch = ___, bty = "n")