AEM 6850

Empirical Methods for Applied Economists

Prof. Ariel Ortiz-Bobea

6 · How to code

Thursday, September 10, 2026

Objectives

By the end of this session you can:

Lay out a project: folders for code, data and output, one script per task, and a master script that runs them in order.
Structure a script: sections, names, comments that say why, named constants, no duplicated code.
Write a function and call it with its own arguments.
Run one function over every file in a folder with lapply().
Stop a script with stopifnot() when a number is not what you wrote down.
Localize a mistake by halving the script, and write the diagnosis down: line, mechanism, fix, proof.

Plan for today

Motivation and philosophy
Organizing projects, scripts, and code
Guiding principles
Writing functions
Checking a script as you write it
Finding a mistake when the number is wrong

Disclaimers

I’m not a computer scientist.
The ideas here reflect my experience and that of others.
Although our work is in R, many of these ideas apply to other languages.

Motivation

Good code…

Increases confidence in your work
Increases the chances other researchers will reproduce, replicate and cite your work
Saves you (and others) time
Reduces the chances of making mistakes
Is easy to understand
May be portable

Coding philosophy

You are writing code for others. That includes yourself at \(t+1\).
A script is written once and read many times.
What you knew when you wrote it is not in the file. Only what you wrote down is.

Two kinds of efficiency

Fast for the machine: fewer operations, less memory, shorter runtime.
Fast for you: you can read it, change it, and find the mistake.
There can be a tradeoff, but in our context writing code that is efficient for you to read and understand dominates.

Organizing a project

A project comes with general instructions
- It could be a readme.txt file
- It should specify the authors, what the project is about, and describe the files
Folders for different types of files
- Don’t mix raw data with “cleaned” data
No single best way of organizing, but this is a good start:
- Folder for code (no subfolders)
- Folder(s) for input(s): raw data you never modify, and cleaned data your scripts create
- Folder(s) for output(s): figures, tables
- Folder for packages (optional)

Example of folder structure

code — all script files, including the RStudio project
data
- raw — data as downloaded from the web, can be large
- clean — data needed for the analysis
output — what shows up in the paper
- figures
- tables
paper
presentations
Others: meeting notes, documentation

Organizing your script files

A script file should accomplish a specific and well-identified task
- Cleaning data, regression analysis, and so on
Never have all your code in a single script file, unless it is a tiny project — your homework, for instance
Script files should have intuitive names and follow a workflow:
- 1.1_clean_economic_data.R
- 1.2_clean_pollution_data.R
- 1.3_prepare_regression_data.R
- 2.1_regression_analysis.R
- 2.2_post_estimation_analysis.R

Good practice to create special script files

A master script file, master.R, that runs all the project’s scripts in the correct order
- Installs and loads all necessary packages
- Simplifies life for whoever reproduces your work
- Typically written at the very end of the project
A file containing user-defined functions, functions.R
- Improves automation when the same function is used in several scripts
- Invites you to write more portable functions for future projects

Code within a script file

# air.R -- days above the PM2.5 standard at one monitor, 2025
# Your Name, NetID

# 0) Setup ----
STANDARD <- 35                        # US EPA 24-hour standard, ug/m3
pm <- read.csv("data/epa_pm25_compton_2025.csv",
               colClasses = c("Site.ID" = "character"))

# 1) One instrument, one month ----

# 2) The numbers ----

# 3) Write results.csv ----

Title block, setup, then one numbered section per task.
# Name ---- is a section. RStudio’s outline lists them.
A script does one job. Two jobs, two scripts.

Two versions of the same script

Both count the days above the standard at one monitor. Both are correct.

STANDARD <- 35                    # US EPA 24-hour standard, ug/m3

pm <- read.csv("data/epa_pm25_compton_2025.csv")
names(pm)[names(pm) == "Daily.Mean.PM2.5.Concentration"] <- "pm25"
pm <- pm[pm$POC == 1, ]           # one instrument

sum(pm$pm25 > STANDARD)

compose <- function(...) { fs <- list(...); function(x) Reduce(\(a, f) f(a), fs, x) }

run <- compose(read.csv, \(d) subset(d, POC == 1),
               \(d) d$Daily.Mean.PM2.5.Concentration, \(v) sum(v > 35))

run("data/epa_pm25_compton_2025.csv")

The second is shorter and more reusable. It is also unreadable at a glance, and the 35 has disappeared into an argument.
Short and clever are not the same thing as clear.

Guiding principles 1/2

Think constantly about reproducibility
- Your script files should run entirely without errors — start with your own computer
- Rely on relative paths, no hard coding
- RStudio projects help. Use them.
Do not repeat yourself
- Do not copy-paste code within a script file. Write a function and use loops.
- Do not copy-paste code across script files. Put the function in its own file and load it when needed.
- Copy-pasting code is a red flag.
Give objects and folders simple and intuitive names
- Good names are short and informative
- Limit upper case and odd characters: MyModified.GDP.var or log_gdp?

Guiding principles 2/2

Annotate your code concisely but generously
- Better to err on the side of too many annotations
- They should help a reader understand the purpose of any significant piece of your code
- Annotate single lines when the code is relatively complex
Make your code look good, which is to say easy to read
- Keep consistent indentation and spacing
- Exploit RStudio’s ability to fold your code
- Avoid deep nesting when possible
- Let your code breathe — leave some blank space
- Favor several lines with short expressions over one long line with several expressions

General tips on the process of writing code

Before starting, visualize the workflow
- “If you don’t know where you’re going, any road’ll take you there” — after Alice in Wonderland, Lewis Carroll
- Have a clear idea of the repetitive tasks
- Think about how the project and the code should be structured when it is done
Then start writing your code interactively
- Working in the console and the editor together
- But do not lose sight of the forest for the trees
Do not run a big task without testing it on a small dataset
Constantly re-check that your code runs entirely, not just in pieces
Leave breadcrumbs so you can find a mistake later

Getting the data

Two files: one monitor in Compton, and every PM2.5 monitor in Los Angeles County, both for 2025. This chunk fetches them if you do not already have them, and does nothing if you do.

dir.create("data", showWarnings = FALSE)

f <- "data/epa_pm25_compton_2025.csv"
if (!file.exists(f)) {
  download.file(paste0("https://arielortizbobea.github.io/aem6850/",
                       "fall-2026/sessions/", f), f)
}

f <- "data/epa_pm25_la_county_2025.csv"
if (!file.exists(f)) {
  download.file(paste0("https://arielortizbobea.github.io/aem6850/",
                       "fall-2026/sessions/", f), f)
}

list.files("data", pattern = "^epa_.*\\.csv$")

#> [1] "epa_pm25_compton_2025.csv"   "epa_pm25_la_county_2025.csv"

dir.create(), then download.file() only if the file is missing.
A script that fetches its own inputs runs on a machine that has never seen the data.

Building the Compton dataframe

Everything below runs off the Compton file, read here so this page and its script stand on their own.

pm <- read.csv("data/epa_pm25_compton_2025.csv",
               colClasses = c("Site.ID" = "character"))
names(pm)[names(pm) == "Daily.Mean.PM2.5.Concentration"] <- "pm25"
names(pm)[names(pm) == "Local.Site.Name"]                <- "site"
pm$date <- as.Date(pm$Date, format = "%m/%d/%Y")
pm <- pm[pm$POC == 1 & pm$AQS.Parameter.Code == 88101, ]

nrow(pm)   # one regulatory instrument, the days it ran

#> [1] 301

Read with the identifier protected, rename, convert the date, keep one instrument. 301 rows.
Every EPA file needs these five lines. A function lets you type them once.

Naming things and named constants

STANDARD <- 35   # US EPA 24-hour PM2.5 standard, ug/m3

sum(pm$pm25 > STANDARD)   # days above the standard, 2025

#> [1] 10

Short, lower case, underscores. Verbs for functions. Never T for TRUE.
A number that means something gets a name at the top, in capitals.
Comments say why. The code already says what.

Write code that is easy to check

Small pieces, one job each.
Intermediate values in named objects you can print.
No step that changes three things at once.
Code you cannot check is code you cannot defend.

`function()`: naming a block of code

days_above <- function(x, standard = 35) {
  sum(x > standard)
}

dec <- pm[format(pm$date, "%m") == "12", ]
days_above(dec$pm25)                  # December, at the 35 standard

#> [1] 5

days_above(dec$pm25, standard = 50)   # the same days, a higher bar

#> [1] 3

days_above(c(1, 50))                  # a case you can check by eye

#> [1] 1

function(x, standard = 35): the arguments. = 35 is a default.
The value is the last expression. return() is for leaving early.
December: 5 days above 35, 3 above 50. c(1, 50): 1, as it must be.

Two versions of the same function

Both return the number of days above the standard.

days_above <- function(x, standard = 35) {
  sum(x > standard)
}

da <- function(v, s = 35, na = TRUE, f = NULL) {
  if (!is.null(f)) v <- f(v)
  length(which(if (na) v[!is.na(v)] > s else v > s))
}

The name does not say what it returns. Four arguments where one is needed, three of them hidden at the call.
A function is easy to read when its name says what it returns and its body fits in your head.

Local names and global names

count_pos <- function(x) {
  n <- sum(x > 0)   # n exists inside the call
  n
}
count_pos(c(-1, 2, 3))

#> [1] 2

n                     # and not outside it

#> Error: object 'n' not found

Names made inside a function die at the closing brace.
So a function can use short names without touching the workspace.

A function can read a global name

standard <- 35
share_above <- function(x) sum(x > standard) / length(x)
share_above(dec$pm25)

#> [1] 0.1612903

rm(standard)
share_above(dec$pm25)   # worked until the global was gone

#> Error in share_above(dec$pm25): object 'standard' not found

A function can read a global name. It works until the global is gone.
Everything a function needs comes in through its arguments.

Writing the five lines once

read_epa <- function(path) {
  d <- read.csv(path, colClasses = c("Site.ID" = "character"))
  names(d)[names(d) == "Daily.Mean.PM2.5.Concentration"] <- "pm25"
  names(d)[names(d) == "Local.Site.Name"]                <- "site"
  d$date <- as.Date(d$Date, format = "%m/%d/%Y")
  d[d$POC == 1 & d$AQS.Parameter.Code == 88101, ]
}

files  <- list.files("data", pattern = "^epa_.*\\.csv$", full.names = TRUE)
frames <- lapply(files, read_epa)
names(frames) <- basename(files)

la      <- frames[["epa_pm25_la_county_2025.csv"]]
la_jan  <- la[format(la$date, "%m") == "01", ]
by_site <- split(la_jan, la_jan$site)

sapply(by_site, nrow)

#>                        Compton        Lancaster - Fairgrounds 
#>                             31                             31 
#> Long Beach-Route 710 Near Road  Los Angeles-North Main Street 
#>                             26                             30 
#>                       Pasadena                 Pico Rivera #2 
#>                             11                             11 
#>                         Reseda             Signal Hill (LBSH) 
#>                             11                              6

The five lines, once, with path as the argument. A change happens in one place.
lapply() over every EPA file, split() into one frame per site.
Eight sites, not eleven: three run no regulatory instrument.

# 1. Write mean_above(x, standard = 35): the mean of the readings in x
#    that are above the standard.
# 2. Run it on every site in by_site with sapply().
# 3. One number to compare with the room: Pasadena.
# 4. One site returns NaN. Say why in one sentence, and whether that is
#    wrong.

mean_above <- function(x, standard = 35) mean(x[x > standard])
round(sapply(by_site, function(s) mean_above(s$pm25)), 1)

#>                        Compton        Lancaster - Fairgrounds 
#>                           46.9                            NaN 
#> Long Beach-Route 710 Near Road  Los Angeles-North Main Street 
#>                           39.2                           82.1 
#>                       Pasadena                 Pico Rivera #2 
#>                           58.1                           40.1 
#>                         Reseda             Signal Hill (LBSH) 
#>                           45.1                           49.3

Pasadena 58.1. Lancaster - Fairgrounds NaN: no reading above 35, so the mean of nothing. Honest.

`stopifnot()`: a check that stops the script

jan <- pm[format(pm$date, "%m") == "01", ]
stopifnot(nrow(pm) == 301)
stopifnot(nrow(jan) == 31)

Silent when every condition holds. Stops the script when one fails.
Several conditions in one call. Right after the step it guards.
The message names the condition that failed.

Test small, then run big

days_above(c(1, 50))                              # a case you can check by eye

#> [1] 1

nrow(read_epa(files[1]))                          # one file before the folder

#> [1] 301

system.time(lapply(files, read_epa))              # the whole folder, timed

#>    user  system elapsed 
#>   0.042   0.002   0.044

Never run the big job untested. A function gets a case you can check by eye.
One file before the folder. system.time() around the big call.

Restart R and source the script

# Session > Restart R   (Ctrl+Shift+F10 on Windows, Cmd+Shift+0 on a Mac)
ls()                    # character(0), or it was not a clean restart
getwd()                 # where "data/..." is looked for
list.files("data")      # what is there
source("air.R")         # your script, top to bottom, in a clean session

Restart R, check ls() is empty, then source(). rm() is not a restart.
Paths are relative to getwd(). The .Rproj file sets it when opened.
“File does not exist”: check getwd() and list.files("data") first.

Debugging as a strategy

Read the error before changing anything.
Reproduce it on the smallest input that still fails.
One hypothesis at a time. Change one thing per run.
Then the tool: traceback(), print(), halving.

`traceback()`: where it died

frames <- lapply(c("data/epa_pm25_compton_2025.csv", "data/nope.csv"), read_epa)
traceback()

Error in file(file, "rt") : cannot open the connection
In addition: Warning message:
In file(file, "rt") :
  cannot open file 'data/nope.csv': No such file or directory

5: file(file, "rt")
4: read.table(file = file, header = header, sep = sep, quote = quote, 
       dec = dec, fill = fill, comment.char = comment.char, ...)
3: read.csv(path, colClasses = c(Site.ID = "character"))
2: FUN(X[[i]], ...)
1: lapply(c("data/epa_pm25_compton_2025.csv", "data/nope.csv"), 
       read_epa)

Read it from the bottom: your call, your function, then R’s.
The error is in the line you wrote at 1. The warning names the file.

`print()` and `cat()`: how a value evolved

for (site in names(by_site)[1:3]) {
  s <- by_site[[site]]
  cat(site, ": ", nrow(s), " days, ", days_above(s$pm25), " above\n", sep = "")
}

#> Compton: 31 days, 4 above
#> Lancaster - Fairgrounds: 31 days, 0 above
#> Long Beach-Route 710 Near Road: 26 days, 3 above

print(c(days = nrow(jan), above = days_above(jan$pm25)))

#>  days above 
#>    31     4

cat() pastes text and numbers into one line. print() shows an object.
A named vector prints label and value together.
Take the line out once the question is answered.

Halving: a wrong number that never errored

No error, no warning, and the number is wrong. Nothing to trace.
Write down what each step should produce: rows, dates, a count.
Comment out half. Run. Compare. Halve the half that disagrees.
Repeat until the mistake has nowhere to hide. Usually one line.

Documenting a diagnosis

# Line:      share_above <- function(x) sum(x > standard) / length(x)
# Mechanism: standard is not an argument; the function reads it from the
#            workspace and fails in a fresh session
# Fix:       the standard comes in as an argument, with a default
share_above <- function(x, standard = 35) sum(x > standard) / length(x)
# Proof:
stopifnot(share_above(c(1, 50)) == 0.5)

Line. Mechanism. Fix. Proof, as a check that runs.
Cannot fix it? The same four lines, fix left blank. Never an edited check.

Practice

Eight exercises with solutions, on the session page and in the live script.

Quick reference

On the page: every command on this page in one table, and the habits next to it.