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
library(tidyr)
library(tidycensus)
library(readr)        # or readxl if your file is Excel
library(ggplot2)
library(tigris)
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(biscale)
library(cowplot)

The Data Base

For this project, I’m using the NCES (National Center for Education Statistics) school-level dataset. I use this database regularly at work, and it really is an incredible resource. You can customize the exact variables you want, download only the years you need, and pull a mix of school demographics, SES indicators, and staffing information.

This flexibility makes it easy to work with, especially when you’re connecting school-level characteristics to community-level indicators. I work for a teen crisis hotline and I use this dataset to see how many of the public Arizona schools we are connected with.

For this tutorial, I focus on the relationship between pupil–teacher ratio and median household income. I chose these two variables because they speak to both school-level resources and the broader socioeconomic context of the communities those schools serve. I was curious whether lower-income areas might also be associated with higher pupil–teacher ratios, which could signal differences in staffing capacity.

Step 1: Preparing the Dataset.

NCES files contain a lot of information, so the first step is deciding which fields are relevant. Every data set will automatically come with the school’s name, so to add to that we need:

-Latitude and longitude (required to determine census tracts) -Pupil–teacher ratio -Any other school attributes you may want later

Coordinate information is essential because census tract assignment relies on geographic location, not the school name. I also added the addresses to make sure the latitude and longitude worked.

Step 2: Cleaning the dataset.

In an excel file I filtered out:

-Closed and inactive schools. -Vocational schools. -Schools with 0 enrolled students.

Step 3: Loading and Cleaning Column Names

I totally forgot to clean the column names in the excel file, so I cleaned it after I added the dataset.

nces <- read.csv("nces_pupil_data - ELSI Export.csv")

Here I check the column names so I can correctly copy and paste them while changing it.

glimpse(nces)
## Rows: 2,064
## Columns: 14
## $ School.Name                                                      <chr> "A J …
## $ State.Name..Public.School..Latest.available.year                 <chr> "Ariz…
## $ Location.Address.1..Public.School..2023.24                       <chr> "855 …
## $ Location.City..Public.School..2023.24                            <chr> "NOGA…
## $ Location.State.Abbr..Public.School..2023.24                      <chr> "AZ "…
## $ Location.ZIP..Public.School..2023.24                             <int> 85621…
## $ Pupil.Teacher.Ratio..Public.School..2023.24                      <chr> "19.1…
## $ Full.Time.Equivalent..FTE..Teachers..Public.School..2023.24      <chr> "16",…
## $ Total.Students..All.Grades..Excludes.AE...Public.School..2023.24 <chr> "306"…
## $ Start.of.Year.Status..Public.School..2023.24                     <chr> "1-Op…
## $ Updated.Status..Public.School..2023.24                           <chr> "1-Op…
## $ Shared.Time.School..Public.School..2023.24                       <chr> "2-No…
## $ Latitude..Public.School..2023.24                                 <dbl> 31.34…
## $ Longitude..Public.School..2023.24                                <dbl> -110.…

I needed to change latitude and longitude specifically to match the format functions like st_as_sf() expect.

nces <- nces %>%
  rename(
    lat    = Latitude..Public.School..2023.24, 
    long    = Longitude..Public.School..2023.24 
  )

Assigning Census Tracts to Each School

Step 1: Convert NCES Data Into a Spatial Object

Here, I turn the NCES dataset into an sf object, which basically means R now recognizes each school as a geographic point based on its latitude and longitude.

-“Coords” tells R which columns represent point coordinates -“crs = 4326” sets the coordinate reference system (standard for GPS data) -“remove = FALSE” keeps the original lat/long columns in case I need them later

options(tigris_use_cache = TRUE)

nces_geo <- st_as_sf(
  nces,
  coords = c("long", "lat"),
  crs = 4326,                  
  remove = FALSE               
)

Step 2: Download Arizona Census Tract Shapefiles

This pulls the official census tract boundaries for Arizona from the Census Bureau. Each tract polygon includes a unique GEOID, which is what we’ll attach to each school.

tracts_az <- tracts(
  state = "AZ",
  year  = 2022,
  cb    = TRUE
)

Step 3: Matching Coordinate Systems

In spatial analysis, all layers need to be in the same coordinate system before you join them.

This line transforms the school point data to match the CRS of the census tract polygons. That way, R can correctly figure out which point falls inside which tract.

nces_points <- st_transform(nces_geo, st_crs(tracts_az ))

Step 4: Joining Schools to Census Tracts

For each school, R checks which census tract polygon it falls inside. It then attaches that tract’s GEOID to the school record

nces_with_geoid <- st_join(
  nces_points,
  tracts_az[, c("GEOID")], 
  left = TRUE
)

Merging ACS Median Income

Now that each school has a census tract assigned, the next step is to bring in census data. For this project, I’m using median household income from the American Community Survey (ACS).

Step 1: Pull ACS income data for Arizona.

Here I use the get_acs() function from tidycensus to download median household income for every census tract in Arizona. Then I used the GEOID column to join the ACS income data to each school.

census_api_key("85e2a7edeb3f23962026aae115c869af1e7c58bf")
## To install your API key for use in future sessions, run this function with `install = TRUE`.
acs_income <- get_acs(
  geography = "tract",
  state     = "AZ",
  year      = 2022,
  survey    = "acs5",
  variables = c(income = "B19013_001"),
  geometry  = FALSE
)
## Getting data from the 2018-2022 5-year ACS
merged <- nces_with_geoid %>%
  left_join(acs_income %>% select(GEOID, income = estimate),
            by = "GEOID")

Step 2: Filter out columns

In the NCES dataset, missing pupil–teacher ratios are coded with a dash “—”. I convert those to NA and then turn the column into a numeric value.

After that, I filter out: -Schools without a valid pupil–teacher ratio. -Schools located in tracts without income data.

I could’ve done this with cleaning the data in the excel file, however I thought maybe I would want to do something with the schools that didn’t have a pupil-teacher ratio but I decided against it.

merged_cleaned <- merged %>%
  mutate(
    ptratio = na_if(Pupil.Teacher.Ratio..Public.School..2023.24, "—"),
    ptratio = as.numeric(ptratio)
  )
## Warning: There was 1 warning in `stopifnot()`.
## ℹ In argument: `ptratio = as.numeric(ptratio)`.
## Caused by warning:
## ! NAs introduced by coercion
merged_cleaned <- merged_cleaned %>%
  filter(
    !is.na(ptratio),
    !is.na(income))

Step 4 - Aggregate to the Census Tract Level

Since each census tract may contain multiple schools, I summarize the data to create one row per tract.

For each tract, I compute: -Mean pupil–teacher ratio across all schools within the tract

tract_data <- merged_cleaned %>%
  group_by(GEOID) %>%
  summarise(
    income_tract  = first(na.omit(income)),
    ptratio_mean  = mean(ptratio, na.rm = TRUE),
    .groups = "drop"
  )

Step 5: Join the aggregated data back to the Arizona tract shapes

Finally, I attach the income and pupil–teacher ratio values back onto the census tract polygons. Dropping geometry from the aggregated table ensures the join works correctly.

After joining, I filter out any tracts missing either income or pupil–teacher ratio so the map will display cleanly.

tract_data_df <- tract_data %>% st_drop_geometry()

tract_income_ptr_sf <- tracts_az %>%
  left_join(tract_data_df, by = "GEOID") %>%
  filter(!is.na(income_tract),
         !is.na(ptratio_mean))

Creating a Bivariate Map

I am also in PAF 515 at the moment and it’s been helpful taking both classes and seeing different things you can do with census data. In that class we used a bivariate choropleth map and I absolutely loved how it showed information of how two variables interact.

Step 1 - Create bivariate classes

The bi_class() function from the biscale package creates combined categories based on both variables:

-“x” is median income. -“y” is mean pupil–teacher ratio. -“style =”quantile”” splits each variable into equally sized groups -“dim = 3” creates a 3×3 grid of combinations (low/medium/high for each variable)

The result is a new column called bi_class that encodes which income/ratio combination each tract falls into (for example, low income + high ratio vs. high income + low ratio).

tract_income_ptr_sf <- bi_class(
  tract_income_ptr_sf,
  x     = income_tract,
  y     = ptratio_mean,
  style = "quantile",
  dim   = 3
)

Step 2: Draw the Bivariate Choropleth Map

Here I use geom_sf() to plot each census tract polygon and fill it according to its bi_class value.

This map lets me see where low-income tracts with high pupil–teacher ratios cluster, and how that pattern compares across the state.

bivar_map <- ggplot() +
  geom_sf(
    data  = tract_income_ptr_sf,
    aes(geometry = geometry, fill = bi_class),
    color        = "white",
    size         = 0.1,
    show.legend  = FALSE
  ) +
  bi_scale_fill(pal = "GrPink", dim = 3) +
  labs(
    title    = "Median Income and Mean Pupil–Teacher Ratio by Census Tract",
    subtitle = "Arizona (tracts with NCES school data)"
  ) +
  bi_theme(base_size = 10)
bivar_map

Step 3: Adding a Legend

Each color block shows how a particular combination of income and ratio is encoded visually.

-The x-axis of the legend represents income (low to high) -The y-axis represents pupil–teacher ratio (low to high)

bivar_legend <- bi_legend(
  pal  = "GrPink",
  dim  = 3,
  xlab = "Median Income (low → high)",
  ylab = "Mean Pupil–Teacher Ratio (low → high)",
  size = 8
)

Step 4: Combining the Map and Legend

Finally, I use cowplot to combine the main map and the legend into a single figure. I position the legend in the lower corner so it’s easy to read without covering important parts of the map.

bivar_final <- ggdraw() +
  draw_plot(bivar_map) +
  draw_plot(
    bivar_legend,
    x      = 0.05,  # adjust position as needed
    y      = 0.02,
    width  = 0.25,
    height = 0.25
  )

bivar_final

# Let’s Talk About The Findings! Looking at the map, one of the first things that stands out is how many more census tracts exist in urban areas like the Phoenix Valley and Tucson. Because of that difference, I realized that if my goal was to get a more granular picture of school staffing patterns and neighborhood income, it probably would have been better to focus on a single county rather than mapping the entire state. That would allow the variation between tracts to show up more clearly instead of being visually diluted by the rural space. It’s something I’ll keep in mind for next time.

Even so, the statewide map still gives us useful information. In the rural areas, I was surprised by the number of tracts that ended up in the grey category, indicating low-income and low-ratio tracts. That actually makes sense, rural areas often have fewer teachers, but they also tend to have fewer students, which can produce a comparatively low pupil–teacher ratio even without large staffing resources. It doesn’t automatically mean those schools are better resourced, it just reflects smaller enrollment.

I was expecting more of the bright blue tracts. These represent areas with low median income and high pupil–teacher ratios, which signals a potential disparity. In lower-income regions, students may be more likely to attend schools where teachers are responsible for a higher number of students on average. Even in a high-level visualization like this, that stands out immediately.