Code
# Run this once if the packages are not already installed:
# install.packages("tidyverse")
library(tidyverse)New York City is made up of neighborhoods with different population sizes, income levels, commuting patterns, and economic conditions. Census-tract data provide a useful way to explore how these characteristics vary across the five boroughs.
This vignette demonstrates how to use the Tidyverse to import, inspect, clean, transform, summarize, and visualize a real-world dataset. The analysis focuses on three questions:
The goal is to demonstrate reproducible data-wrangling and visualization techniques. The analysis is descriptive: it does not establish that one characteristic causes another.
The dataset contains demographic, economic, employment, and commuting information for New York City census tracts. It was originally distributed through Kaggle as New York City Census Data and is accessed here from a public GitHub copy of the CSV file.
The source file contains census-tract identifiers, county and borough names, total population, income measures, poverty measures, employment information, and commuting characteristics. The original example reports 2,167 rows and 36 columns. The code below checks the dimensions again when the file is loaded.
Unit of analysis: one row represents a census tract. Therefore, borough averages in this vignette are averages across tracts, not population-weighted estimates for all borough residents.
After completing this vignette, a reader should be able to:
readr::read_csv();glimpse() and summarise();is.na() and across();select();filter() and drop_na();mutate() and case_when();group_by() and summarise();arrange() and slice_*(); andggplot2.The tidyverse package loads a collection of packages designed to work together for data science. This vignette mainly uses readr, dplyr, tidyr, ggplot2, and forcats.
If the packages are not installed, run the installation command once. It is included as a comment so that knitting the document does not attempt to install packages automatically.
# Run this once if the packages are not already installed:
# install.packages("tidyverse")
library(tidyverse)read_csv() reads a comma-separated file and returns a tibble. A tibble is a modern form of data frame that prints useful information without displaying every row and column.
data_url <- paste0("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-607/main/Tidyverse/nyc_census_tracts.csv")
nyc_raw <- read_csv(data_url, show_col_types = FALSE)
# Check the number of rows and columns.
dim(nyc_raw)#> [1] 2167 36
# Display the column names.
names(nyc_raw)#> [1] "CensusTract" "County" "Borough" "TotalPop"
#> [5] "Men" "Women" "Hispanic" "White"
#> [9] "Black" "Native" "Asian" "Citizen"
#> [13] "Income" "IncomeErr" "IncomePerCap" "IncomePerCapErr"
#> [17] "Poverty" "ChildPoverty" "Professional" "Service"
#> [21] "Office" "Construction" "Production" "Drive"
#> [25] "Carpool" "Transit" "Walk" "OtherTransp"
#> [29] "WorkAtHome" "MeanCommute" "Employed" "PrivateWork"
#> [33] "PublicWork" "SelfEmployed" "FamilyWork" "Unemployment"
# Inspect variable types and example values.
glimpse(nyc_raw)#> Rows: 2,167
#> Columns: 36
#> $ CensusTract <dbl> 36005000100, 36005000200, 36005000400, 36005001600, 36…
#> $ County <chr> "Bronx", "Bronx", "Bronx", "Bronx", "Bronx", "Bronx", …
#> $ Borough <chr> "Bronx", "Bronx", "Bronx", "Bronx", "Bronx", "Bronx", …
#> $ TotalPop <dbl> 7703, 5403, 5915, 5879, 2591, 8516, 4774, 150, 5355, 3…
#> $ Men <dbl> 7133, 2659, 2896, 2558, 1206, 3301, 2130, 109, 2338, 1…
#> $ Women <dbl> 570, 2744, 3019, 3321, 1385, 5215, 2644, 41, 3017, 164…
#> $ Hispanic <dbl> 29.9, 75.8, 62.7, 65.1, 55.4, 61.1, 62.3, 0.0, 76.5, 6…
#> $ White <dbl> 6.1, 2.3, 3.6, 1.6, 9.0, 1.6, 0.2, 52.0, 1.5, 0.0, 1.6…
#> $ Black <dbl> 60.9, 16.0, 30.7, 32.4, 29.0, 31.1, 36.5, 48.0, 18.9, …
#> $ Native <dbl> 0.2, 0.0, 0.0, 0.0, 0.0, 0.3, 1.0, 0.0, 0.0, 0.0, 0.0,…
#> $ Asian <dbl> 1.6, 4.2, 0.3, 0.0, 2.1, 3.3, 0.0, 0.0, 3.0, 0.0, 0.0,…
#> $ Citizen <dbl> 6476, 3639, 4100, 3536, 1557, 5436, 3056, 41, 2509, 14…
#> $ Income <dbl> NA, 72034, 74836, 32312, 37936, 18086, 14479, NA, 1722…
#> $ IncomeErr <dbl> NA, 13991, 8407, 6859, 3771, 3694, 1901, NA, 6097, 522…
#> $ IncomePerCap <dbl> 2440, 22180, 27700, 17526, 17986, 12023, 9781, 49365, …
#> $ IncomePerCapErr <dbl> 373, 2206, 2449, 2945, 2692, 2346, 1460, 2072, 1714, 1…
#> $ Poverty <dbl> NA, 20.0, 13.2, 26.3, 37.1, 53.2, 54.2, 0.0, 50.5, 52.…
#> $ ChildPoverty <dbl> NA, 20.7, 23.6, 35.9, 31.5, 67.7, 68.3, 0.0, 62.4, 64.…
#> $ Professional <dbl> NA, 28.7, 32.2, 19.1, 35.4, 14.5, 20.8, 63.7, 15.9, 19…
#> $ Service <dbl> NA, 17.1, 23.4, 36.1, 20.9, 41.1, 36.0, 36.3, 36.6, 35…
#> $ Office <dbl> NA, 23.9, 24.9, 26.2, 26.2, 16.7, 33.7, 0.0, 24.2, 25.…
#> $ Construction <dbl> NA, 8.0, 9.0, 4.9, 6.6, 7.1, 3.4, 0.0, 8.6, 4.9, 1.7, …
#> $ Production <dbl> NA, 22.3, 10.5, 13.8, 11.0, 20.6, 6.1, 0.0, 14.7, 13.5…
#> $ Drive <dbl> NA, 44.8, 41.3, 37.2, 19.2, 19.6, 5.9, 0.0, 12.6, 14.0…
#> $ Carpool <dbl> NA, 13.7, 10.0, 5.3, 5.3, 7.0, 0.0, 0.0, 2.8, 1.7, 3.0…
#> $ Transit <dbl> NA, 38.6, 44.6, 45.5, 63.9, 68.2, 74.5, 100.0, 62.5, 6…
#> $ Walk <dbl> NA, 2.9, 1.4, 8.6, 3.0, 4.3, 14.0, 0.0, 17.7, 18.0, 7.…
#> $ OtherTransp <dbl> NA, 0.0, 0.5, 1.6, 2.4, 1.0, 1.5, 0.0, 1.8, 0.0, 0.7, …
#> $ WorkAtHome <dbl> NA, 0.0, 2.1, 1.7, 6.2, 0.0, 4.1, 0.0, 2.7, 1.6, 0.5, …
#> $ MeanCommute <dbl> NA, 43.0, 45.0, 38.8, 45.4, 46.0, 42.7, NA, 35.5, 42.8…
#> $ Employed <dbl> 0, 2308, 2675, 2120, 1083, 2508, 1191, 113, 1691, 1102…
#> $ PrivateWork <dbl> NA, 80.8, 71.7, 75.0, 76.8, 71.0, 74.2, 62.8, 85.1, 86…
#> $ PublicWork <dbl> NA, 16.2, 25.3, 21.3, 15.5, 21.3, 16.1, 37.2, 8.3, 8.5…
#> $ SelfEmployed <dbl> NA, 2.9, 2.5, 3.8, 7.7, 7.7, 9.7, 0.0, 6.1, 4.5, 11.0,…
#> $ FamilyWork <dbl> NA, 0.0, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, …
#> $ Unemployment <dbl> NA, 7.7, 9.5, 8.7, 19.2, 17.2, 18.9, 0.0, 9.4, 15.2, 1…
The original dataset contains 2,167 census-tract records and 36 variables. The dim() and glimpse() commands let us verify the structure rather than assume that the downloaded file has not changed.
Missing data can affect summaries and plots. Before cleaning the dataset, it is useful to determine which columns contain missing values.
The code below combines summarise(), across(), and pivot_longer():
summarise() creates a one-row summary;across() applies the same calculation to every column;is.na() identifies missing values; andpivot_longer() changes the summary from a wide layout to a tidy two-column layout.missing_summary <- nyc_raw %>%
summarise(
across(
everything(),
~ sum(is.na(.x))
)
) %>%
pivot_longer(
cols = everything(),
names_to = "variable",
values_to = "missing_count"
) %>%
mutate(
missing_percent = round(
100 * missing_count / nrow(nyc_raw),
2
)
) %>%
filter(missing_count > 0) %>%
arrange(desc(missing_count))
missing_summary#> # A tibble: 28 × 3
#> variable missing_count missing_percent
#> <chr> <int> <dbl>
#> 1 Income 66 3.05
#> 2 IncomeErr 66 3.05
#> 3 MeanCommute 61 2.81
#> 4 ChildPoverty 60 2.77
#> 5 IncomePerCap 46 2.12
#> 6 IncomePerCapErr 46 2.12
#> 7 Professional 43 1.98
#> 8 Service 43 1.98
#> 9 Office 43 1.98
#> 10 Construction 43 1.98
#> # ℹ 18 more rows
This table reports the number and percentage of missing values for each affected variable. It helps identify which fields need attention before analysis.
The analysis does not require all 36 columns. select() keeps only the variables needed for the research questions. filter() removes records without the borough, income, or population information required for the main analysis.
The income and population fields are also checked for positive values. The selected poverty and commute fields are retained even if they contain missing values; they will be handled when the relevant analysis is performed.
nyc <- nyc_raw %>%
select(
CensusTract,
County,
Borough,
TotalPop,
Income,
IncomePerCap,
Poverty,
ChildPoverty,
MeanCommute,
Transit,
Walk,
Drive,
Unemployment
) %>%
filter(
!is.na(Borough),
!is.na(Income),
!is.na(TotalPop),
Income > 0,
TotalPop > 0
) %>%
mutate(
Borough = as.factor(Borough)
)
# Confirm the cleaned data structure.
glimpse(nyc)#> Rows: 2,101
#> Columns: 13
#> $ CensusTract <dbl> 36005000200, 36005000400, 36005001600, 36005001900, 36005…
#> $ County <chr> "Bronx", "Bronx", "Bronx", "Bronx", "Bronx", "Bronx", "Br…
#> $ Borough <fct> Bronx, Bronx, Bronx, Bronx, Bronx, Bronx, Bronx, Bronx, B…
#> $ TotalPop <dbl> 5403, 5915, 5879, 2591, 8516, 4774, 5355, 3016, 4778, 529…
#> $ Income <dbl> 72034, 74836, 32312, 37936, 18086, 14479, 17226, 20153, 1…
#> $ IncomePerCap <dbl> 22180, 27700, 17526, 17986, 12023, 9781, 11493, 10317, 89…
#> $ Poverty <dbl> 20.0, 13.2, 26.3, 37.1, 53.2, 54.2, 50.5, 52.6, 52.7, 8.5…
#> $ ChildPoverty <dbl> 20.7, 23.6, 35.9, 31.5, 67.7, 68.3, 62.4, 64.9, 63.1, 6.5…
#> $ MeanCommute <dbl> 43.0, 45.0, 38.8, 45.4, 46.0, 42.7, 35.5, 42.8, 44.0, 47.…
#> $ Transit <dbl> 38.6, 44.6, 45.5, 63.9, 68.2, 74.5, 62.5, 64.7, 83.6, 48.…
#> $ Walk <dbl> 2.9, 1.4, 8.6, 3.0, 4.3, 14.0, 17.7, 18.0, 7.1, 2.0, 14.6…
#> $ Drive <dbl> 44.8, 41.3, 37.2, 19.2, 19.6, 5.9, 12.6, 14.0, 5.1, 42.8,…
#> $ Unemployment <dbl> 7.7, 9.5, 8.7, 19.2, 17.2, 18.9, 9.4, 15.2, 10.6, 12.8, 6…
# Count the records retained for each borough.
nyc %>%
count(Borough, name = "tract_count") %>%
arrange(desc(tract_count))#> # A tibble: 5 × 2
#> Borough tract_count
#> <fct> <int>
#> 1 Brooklyn 748
#> 2 Queens 642
#> 3 Bronx 329
#> 4 Manhattan 275
#> 5 Staten Island 107
This cleaning step creates an analysis dataset containing records with the minimum information needed for the income and population comparisons. It does not remove a tract merely because a secondary variable, such as commute time, is missing.
A numeric poverty percentage can be easier to interpret when grouped into meaningful categories. mutate() creates a new column, while case_when() assigns a label based on conditions.
The categories below are descriptive groupings selected for this example; they are not official government classifications.
nyc <- nyc %>%
mutate(
PovertyGroup = case_when(
is.na(Poverty) ~ "Not reported",
Poverty < 10 ~ "Below 10%",
Poverty < 20 ~ "10% to below 20%",
TRUE ~ "20% or higher"
),
PovertyGroup = factor(
PovertyGroup,
levels = c(
"Below 10%",
"10% to below 20%",
"20% or higher",
"Not reported"
)
)
)
nyc %>%
count(PovertyGroup, name = "tract_count") %>%
arrange(PovertyGroup)#> # A tibble: 3 × 2
#> PovertyGroup tract_count
#> <fct> <int>
#> 1 Below 10% 571
#> 2 10% to below 20% 675
#> 3 20% or higher 855
The new PovertyGroup variable can be used to compare income across broad poverty categories. Keeping missing poverty values as "Not reported" also makes the missing-data category visible rather than silently treating it as zero.
group_by() divides the data into borough groups. summarise() then calculates statistics within each group. The na.rm = TRUE argument tells the summary functions to ignore missing values for the particular variable being summarized.
The income statistics below describe the distribution of census-tract income values within each borough. They should not be interpreted as population-weighted borough income estimates.
borough_summary <- nyc %>%
group_by(Borough) %>%
summarise(
census_tracts = n(),
mean_tract_income = mean(Income, na.rm = TRUE),
median_tract_income = median(Income, na.rm = TRUE),
mean_population = mean(TotalPop, na.rm = TRUE),
mean_poverty_percent = mean(Poverty, na.rm = TRUE),
mean_commute_minutes = mean(MeanCommute, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(desc(median_tract_income))
borough_summary#> # A tibble: 5 × 7
#> Borough census_tracts mean_tract_income median_tract_income mean_population
#> <fct> <int> <dbl> <dbl> <dbl>
#> 1 Manhattan 275 83352. 81250 5875.
#> 2 Staten Is… 107 72021. 75167 4415.
#> 3 Queens 642 62337. 60629 3583.
#> 4 Brooklyn 748 54265. 48398. 3467.
#> 5 Bronx 329 39311. 33365 4312.
#> # ℹ 2 more variables: mean_poverty_percent <dbl>, mean_commute_minutes <dbl>
The resulting table makes it possible to compare the boroughs using the same set of summary measures. The median is included because it is less sensitive than the mean to unusually high or low tract-level income values.
arrange() sorts rows, and slice_head() and slice_tail() select records from the sorted result. The following code identifies a small number of census tracts at each end of the income distribution.
highest_income_tracts <- nyc %>%
arrange(desc(Income)) %>%
select(
CensusTract,
County,
Borough,
TotalPop,
Income,
Poverty
) %>%
slice_head(n = 10)
lowest_income_tracts <- nyc %>%
arrange(Income) %>%
select(
CensusTract,
County,
Borough,
TotalPop,
Income,
Poverty
) %>%
slice_head(n = 10)
highest_income_tracts#> # A tibble: 10 × 6
#> CensusTract County Borough TotalPop Income Poverty
#> <dbl> <chr> <fct> <dbl> <dbl> <dbl>
#> 1 36061013000 New York Manhattan 3277 244375 7.1
#> 2 36061014200 New York Manhattan 4467 235455 5.9
#> 3 36061016001 New York Manhattan 3861 222222 2
#> 4 36047002100 Kings Brooklyn 4295 208036 7.6
#> 5 36061011401 New York Manhattan 1276 197500 3.2
#> 6 36061011202 New York Manhattan 287 192857 4.5
#> 7 36061006900 New York Manhattan 2604 183750 7.4
#> 8 36061003100 New York Manhattan 2673 183500 13.1
#> 9 36061031703 New York Manhattan 8875 182401 4.3
#> 10 36061012800 New York Manhattan 6226 175259 4.5
lowest_income_tracts#> # A tibble: 10 × 6
#> CensusTract County Borough TotalPop Income Poverty
#> <dbl> <chr> <fct> <dbl> <dbl> <dbl>
#> 1 36047080800 Kings Brooklyn 1256 9829 71.4
#> 2 36047091000 Kings Brooklyn 6094 10746 67.3
#> 3 36047035200 Kings Brooklyn 1035 11011 51.6
#> 4 36061011900 New York Manhattan 1068 11671 55.5
#> 5 36005014300 Bronx Bronx 968 12000 52.2
#> 6 36005016100 Bronx Bronx 4140 12467 51.6
#> 7 36047090800 Kings Brooklyn 5679 12863 62.1
#> 8 36047008500 Kings Brooklyn 6916 14000 52.7
#> 9 36085013301 Richmond Staten Island 1531 14030 68.3
#> 10 36005014701 Bronx Bronx 7023 14334 61
These results show the range of tract-level income values in the cleaned dataset. They are useful for checking whether a small number of unusually high or low observations may influence the overall distribution.
A histogram shows how frequently income values occur across the census tracts. The vertical line marks the median tract-level income in the cleaned dataset.
ggplot() defines the data and aesthetic mapping, while geom_histogram() draws the bars. geom_vline() adds the median reference line.
overall_median_income <- median(nyc$Income, na.rm = TRUE)
ggplot(nyc, aes(x = Income)) +
geom_histogram(
bins = 40,
fill = "steelblue",
color = "white"
) +
geom_vline(
xintercept = overall_median_income,
color = "red",
linewidth = 1,
linetype = "dashed"
) +
labs(
title = "Distribution of Census-Tract Income in New York City",
subtitle = "The dashed line represents the median across included tracts",
x = "Reported income (USD)",
y = "Number of census tracts",
caption = "Source: New York City Census Data; see dataset links above."
) +
theme_minimal()The histogram helps show whether tract-level income values are concentrated in a particular range or whether the distribution has a long tail. The median line provides a reference point for the distribution.
A boxplot displays the median, interquartile range, and potential outlying values for each borough. It is useful for comparing the spread of tract-level income values, rather than comparing only one average.
fct_reorder() orders the borough labels by their median income, making the chart easier to read.
nyc %>%
mutate(
Borough = fct_reorder(Borough, Income, .fun = median, na.rm = TRUE)
) %>%
ggplot(aes(x = Borough, y = Income)) +
geom_boxplot(
fill = "skyblue",
color = "gray25",
outlier.alpha = 0.35
) +
coord_flip() +
scale_y_continuous(labels = scales::label_dollar()) +
labs(
title = "Census-Tract Income Distribution by Borough",
x = "Borough",
y = "Reported income (USD)",
caption = "Each observation represents a census tract."
) +
theme_minimal()The boxplots show both differences in typical income and differences in the amount of variation within each borough. Overlapping boxes indicate that the borough distributions are not completely separated.
A scatterplot can help determine whether census tracts with higher poverty percentages tend to have lower reported income. Each point represents one tract.
geom_smooth(method = "lm") adds a fitted linear trend line. This line summarizes an association in the observed data; it does not establish a causal relationship.
poverty_income_data <- nyc %>%
drop_na(Poverty, Income)
ggplot(
poverty_income_data,
aes(x = Poverty, y = Income)
) +
geom_point(
alpha = 0.45,
color = "steelblue"
) +
geom_smooth(
method = "lm",
se = TRUE,
color = "darkred"
) +
scale_y_continuous(labels = scales::label_dollar()) +
labs(
title = "Census-Tract Poverty and Income",
subtitle = "Each point represents a tract with both measures available",
x = "Population below poverty level (%)",
y = "Reported income (USD)",
caption = "The fitted line describes association, not causation."
) +
theme_minimal()The direction and spread of the points provide a visual indication of whether poverty and income are associated. The width of the scatter also reminds us that poverty alone does not explain all differences in tract-level income.
The dataset also contains MeanCommute, which records mean commute time in minutes. This section uses drop_na() to include only tracts with a reported commute-time value.
commute_summary <- nyc %>%
drop_na(MeanCommute) %>%
group_by(Borough) %>%
summarise(
mean_commute_minutes = mean(MeanCommute),
median_commute_minutes = median(MeanCommute),
tracts_with_commute_data = n(),
.groups = "drop"
) %>%
arrange(desc(mean_commute_minutes))
commute_summary#> # A tibble: 5 × 4
#> Borough mean_commute_minutes median_commute_minutes tracts_with_commute_…¹
#> <fct> <dbl> <dbl> <int>
#> 1 Queens 42.9 42.4 641
#> 2 Bronx 42.6 42.6 327
#> 3 Staten Isl… 42.2 41.9 107
#> 4 Brooklyn 41.9 42 748
#> 5 Manhattan 30.8 30.3 275
#> # ℹ abbreviated name: ¹tracts_with_commute_data
The next plot compares the distribution of tract-level commute times across boroughs.
nyc %>%
drop_na(MeanCommute) %>%
mutate(
Borough = fct_reorder(
Borough,
MeanCommute,
.fun = median
)
) %>%
ggplot(aes(x = Borough, y = MeanCommute)) +
geom_boxplot(
fill = "tan",
color = "gray25"
) +
coord_flip() +
labs(
title = "Mean Commute Time by Borough",
x = "Borough",
y = "Mean commute time (minutes)",
caption = "Each observation represents a census tract with commute data."
) +
theme_minimal()The chart can reveal differences in typical commute time and variation across tracts within each borough. It does not identify the transportation mode or explain why commute times differ.
The pivot_longer() function converts selected summary columns into a long format. This is useful when several measures need to be displayed using the same plotting structure.
The measures in this example have different units, so the plot uses separate panels with independent y-axis scales.
borough_summary_long <- borough_summary %>%
select(
Borough,
median_tract_income,
mean_poverty_percent,
mean_commute_minutes
) %>%
pivot_longer(
cols = -Borough,
names_to = "measure",
values_to = "value"
)
borough_summary_long#> # A tibble: 15 × 3
#> Borough measure value
#> <fct> <chr> <dbl>
#> 1 Manhattan median_tract_income 81250
#> 2 Manhattan mean_poverty_percent 17.6
#> 3 Manhattan mean_commute_minutes 30.8
#> 4 Staten Island median_tract_income 75167
#> 5 Staten Island mean_poverty_percent 13.3
#> 6 Staten Island mean_commute_minutes 42.2
#> 7 Queens median_tract_income 60629
#> 8 Queens mean_poverty_percent 14.0
#> 9 Queens mean_commute_minutes 42.9
#> 10 Brooklyn median_tract_income 48398.
#> 11 Brooklyn mean_poverty_percent 21.6
#> 12 Brooklyn mean_commute_minutes 41.9
#> 13 Bronx median_tract_income 33365
#> 14 Bronx mean_poverty_percent 29.4
#> 15 Bronx mean_commute_minutes 42.6
ggplot(
borough_summary_long,
aes(
x = fct_reorder(Borough, value),
y = value
)
) +
geom_col(fill = "cadetblue") +
coord_flip() +
facet_wrap(
~ measure,
scales = "free_x"
) +
labs(
title = "Selected Census-Tract Measures by Borough",
x = "Borough",
y = "Value"
) +
theme_minimal()This example demonstrates how pivot_longer() makes a wide summary table suitable for faceted plotting. Because income, poverty, and commute time use different units, each panel has its own scale and should be interpreted separately.
Use the summary tables and plots generated above to describe the patterns in the downloaded data. The following points are intended as an interpretation guide; the actual values should be read from the knitted output.
These findings are descriptive and should not be interpreted as causal. The dataset represents census tracts, not individual households, and the borough-level summaries are not population-weighted estimates. Other variables and data from different years could produce different patterns.
This vignette demonstrated a practical Tidyverse workflow using New York City census-tract data. The analysis imported a CSV file, inspected its structure, measured missingness, selected relevant columns, created a poverty category, grouped observations by borough, calculated descriptive summaries, and created visualizations.
The examples highlighted several useful functions:
read_csv()select(), filter(), mutate(), case_when(), group_by(), summarise(), arrange(), count(), drop_na(), and slice_head()pivot_longer()fct_reorder()ggplot(), aes(), geom_histogram(), geom_vline(), geom_boxplot(), geom_point(), geom_smooth(), geom_col(), facet_wrap(), and coord_flip()Together, these functions make it possible to create a readable and reproducible data-analysis workflow. The same approach can be adapted to other public datasets and research questions.