AEM 6850

Empirical Methods for Applied Economists

Prof. Ariel Ortiz-Bobea

Wrangling with dplyr

Optional · not taught in class

Cornell University

What this page is

This page is not on the calendar. It is written, runnable, and about one class meeting of material.

Objectives

By the end of this page you can:

  • Translate the base-R subsetting you already know into dplyr verbs.
  • Build grouped summaries with group_by() and summarize().
  • Read a pipe chain aloud as a sentence.
  • Predict row and group counts before running.
  • Join two tables on a shared key with left_join(), and check the row count before and after.

The first package

A package is R code somebody else wrote, tested, and shared. dplyr is the standard one for wrangling data frames. In the console:

install.packages("dplyr")
  • A minute of console messages is normal. It ends with one line saying the package was installed.
  • Once per machine, like installing R itself. Never in a script.

Load it, every session

Installing put dplyr on the machine. library() puts it into the current R session:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
  • The masked-objects block is information, not an error. dplyr’s filter() and lag() now stand in front of the base ones.
  • Install once per machine. library() at the top of every script.

Rebuild the frame

Reading stays base R: the same read, the same repairs, the same checks.

pm <- read.csv("data/epa_pm25_compton_2025.csv",
               colClasses = c("Site.ID"          = "character",
                              "State.FIPS.Code"  = "character",
                              "County.FIPS.Code" = "character"))
names(pm)[names(pm) == "Daily.Mean.PM2.5.Concentration"] <- "pm25"
pm$date <- as.Date(pm$Date, format = "%m/%d/%Y")

nrow(pm)   # the whole year, both instruments: 662 rows
#> [1] 662

filter(): keep rows

Base R cuts the winter window with brackets. Same subset, first verb:

# base R:  pm[pm$date <= as.Date("2025-02-28") &
#                pm$POC == 1 & pm$AQS.Parameter.Code == 88101, ]

one <- filter(pm,
              date <= as.Date("2025-02-28"),
              POC == 1,
              AQS.Parameter.Code == 88101)

nrow(one)   # must be 59: one row per calendar day
#> [1] 59
  • Data frame first, then conditions. Bare column names: no pm$.
  • Commas mean and. | still spells or, inside one condition.
  • Same rows out: 59, one per day.

select(): keep columns

Base needed two coordinates and quoted names. select() takes bare names, in the order you want them:

# base R:  one[, c("date", "pm25")]

head(select(one, date, pm25), 3)
#>         date pm25
#> 1 2025-01-01 53.2
#> 2 2025-01-02 33.6
#> 3 2025-01-03 27.7

mutate(): add columns

# base R:  one$above35 <- one$pm25 > 35

one <- mutate(one,
              above35 = pm25 > 35,
              month   = format(date, "%b"))

head(select(one, date, pm25, above35, month), 3)
#>         date pm25 above35 month
#> 1 2025-01-01 53.2    TRUE   Jan
#> 2 2025-01-02 33.6   FALSE   Jan
#> 3 2025-01-03 27.7   FALSE   Jan
sum(one$above35)   # must be 4
#> [1] 4
  • name = expression, computed inside the frame: bare column names again.
  • "%b" is the short month name: format strings, in reverse.

arrange(): sort rows

arrange() sorts rows by a column. desc() flips it to descending:

# base R:  head(one[order(-one$pm25), c("date", "pm25")], 5)

worst <- arrange(one, desc(pm25))
head(select(worst, date, pm25), 3)
#>         date pm25
#> 1 2025-01-01 53.2
#> 2 2025-01-09 47.7
#> 3 2025-01-08 44.6
  • January 1 on top at 53.2, the same answer base sorting gave.
  • The cost so far: throwaway objects like worst, or nested calls you read inside out. The pipe fixes this.

The pipe

|> passes whatever is on its left to the function on its right, as that function’s first argument:

nrow(filter(pm, POC == 1))         # nested: read inside out
#> [1] 301
pm |> filter(POC == 1) |> nrow()   # piped: read left to right
#> [1] 301
  • Same computation, same 301. The piped line reads in the order things happen.
  • Read |> as “then”. Verbs take the data frame first so chains work.
  • Older code writes %>%. Read it the same way; write |>.

The chain, end to end

# base R:  head(one[order(-one$pm25), c("date", "pm25")], 5)

pm |>
  filter(date <= as.Date("2025-02-28"), POC == 1, AQS.Parameter.Code == 88101) |>
  arrange(desc(pm25)) |>
  select(date, pm25) |>
  head(5)
#>         date pm25
#> 1 2025-01-01 53.2
#> 2 2025-01-09 47.7
#> 3 2025-01-08 44.6
#> 4 2025-01-10 42.2
#> 5 2025-01-02 33.6
  • Aloud: take the year, keep the winter window, sort worst first, keep two columns, show five.
  • Same five rows as the base sort: 53.2, 47.7, 44.6, 42.2, 33.6.
  • The whole thing reads as one sentence, one verb per line.

group_by() and summarize()

Write down what this should print. Evidence: no gaps in the window, so 31 and 28 days. January has all four event days, so its mean should clear its median with room.

one |>
  group_by(month) |>
  summarize(days        = n(),
            mean_pm25   = mean(pm25),
            median_pm25 = median(pm25))
#> # A tibble: 2 x 4
#>   month  days mean_pm25 median_pm25
#>   <chr> <int>     <dbl>       <dbl>
#> 1 Feb      28      13.4        12.6
#> 2 Jan      31      21.3        16.9

Reading the grouped result

  • Days 31 and 28. n() counts the rows of its group; it works only inside summarize().
  • Medians 16.9 and 12.6. January’s mean sits at 21.3, pulled up by the four event days.
  • “A tibble”: dplyr’s flavor of data frame. Same behavior, more careful printing.
  • Nobody predicted this part: February printed first. Text sorts alphabetically.

Factors put months in order

month.abb   # a vector R ships with
#>  [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"
one <- mutate(one, month = factor(month, levels = month.abb))

one |>
  group_by(month) |>
  summarize(days = n(), median_pm25 = median(pm25))
#> # A tibble: 2 x 3
#>   month  days median_pm25
#>   <fct> <int>       <dbl>
#> 1 Jan      31        16.9
#> 2 Feb      28        12.6
  • The level table carries the order you meant. Every grouped result now follows it.
  • This is the factor’s job: categories with a real order, recorded once.

count(): the counting shorthand

# base R:  table(pm$POC, pm$AQS.Parameter.Code)

pm |> count(POC, AQS.Parameter.Code)
#>   POC AQS.Parameter.Code   n
#> 1   1              88101 301
#> 2   3              88101   7
#> 3   3              88502 354
  • count(x) = group_by(x) + summarize(n = n()).
  • Same three cells as table(): 301, 7, 354. Sum: 662, the whole file.
  • The difference: the answer is a data frame, so it can go through more verbs.

The missing spring, measured

The regulatory series has 301 rows where a full year has 365. A chain you can now read says which months the missing days fall in:

Write down the twelve counts this should print. The evidence: the count() above says the regulatory series has 301 rows, and a full month has 28 to 31 days, so 64 days are missing somewhere. Your twelve numbers must sum to 301. Decide where you think the gap falls before you look.

year <- pm |>
  filter(POC == 1, AQS.Parameter.Code == 88101) |>
  mutate(month = factor(format(date, "%b"), levels = month.abb))

The twelve counts

year |> count(month)
#>    month  n
#> 1    Jan 31
#> 2    Feb 28
#> 3    Mar 30
#> 4    Apr  2
#> 5    May  1
#> 6    Jun 28
#> 7    Jul 31
#> 8    Aug 31
#> 9    Sep 27
#> 10   Oct 31
#> 11   Nov 30
#> 12   Dec 31
  • April 2 and May 1: the sampler was down for the spring. The counts sum to 301, so 64 days are missing from a 365-day year.
  • Smaller gaps: March 30 of 31, June 28 of 30, September 27 of 30.

Reading the coverage table

  • 58 of the 64 missing days are the spring outage. The other six are single days in March, June, and September.
  • The cost, quantified on the page: the missing months were the cleanest stretch of the year by the site’s other instrument. An annual mean over the 301 days leans high.

Reading the weather file

Readings in one table, the weather in another. The weather file carries three lines of station metadata before its header, so skip = 3:

wx <- read.csv("data/la-weather-dec2024-feb2025.csv", skip = 3)
names(wx)[names(wx) == "precipitation_sum..mm."] <- "rain"
wx$date <- as.Date(wx$time)

nrow(one); nrow(wx)   # count BEFORE: 59 days of readings, 90 days of weather
#> [1] 59
#> [1] 90
  • skip = 3 steps over the station metadata. Rename, convert the date.
  • 59 days of readings, 90 days of weather, one date column in each.

left_join(): joining two tables

# base R:  merge(one, wx, by = "date", all.x = TRUE)

both <- left_join(select(one, date, pm25),
                  select(wx, date, rain),
                  by = "date")
nrow(both)            # and after: 59
#> [1] 59
sum(is.na(both$rain))
#> [1] 0
  • by names the shared column. The left table’s rows, both tables’ columns.
  • 59 in, 59 out, no NA: every day found its weather.
  • Write down the row count you expect, then compare. Every time.

Rain days against dry days

The join is what makes the next question one chain:

both |>
  mutate(wet = rain > 0) |>
  group_by(wet) |>
  summarize(days = n(), mean_pm25 = mean(pm25))
#> # A tibble: 2 x 3
#>   wet    days mean_pm25
#>   <lgl> <int>     <dbl>
#> 1 FALSE    48     20.1 
#> 2 TRUE     11      6.43
  • 11 rain days averaged 6.4. The 48 dry days averaged 20.1.
  • The biggest storm, 37.9 mm on February 13, sits under a reading of 3.2.

How a join drops and multiplies rows

partial <- filter(wx, date != as.Date("2025-01-07"))   # pretend one day is missing

nrow(inner_join(one, partial, by = "date"))   # dropped
#> [1] 58
nrow(left_join(one, partial, by = "date"))    # kept, with NA
#> [1] 59
  • inner_join() keeps only keys in both tables: rows vanish, silently.
  • left_join() keeps the left table and fills with NA.
  • A duplicated key on the right multiplies rows instead.
  • Mismatched key types stop with an error. The only failure that does.

The one-line stopwatch

One more habit. Wrap any line in system.time() and R reports how long it took:

system.time(read.csv("data/epa_pm25_compton_2025.csv"))
#>    user  system elapsed 
#>   0.004   0.000   0.004
  • Read elapsed: the seconds you waited. Here, effectively zero.
  • The habit: when a line feels slow, wrap it once and note the number.

Quick reference

On the page: the verb table with base translations, the pipe and its shortcut, the factor recipe, and the things that bite.