# tidyverse gives me readr for import, dplyr for cleaning, and ggplot2 for plots.
# ggrepel places the country labels on the final plot so they do not overlap each other.
library(tidyverse)
library(ggrepel)What Makes a Better Cup? Growing Altitude, Origin, and Coffee Quality
Introduction
Coffee professionals talk constantly about “high-grown” coffee as though altitude were a guarantee of quality. This project tests that claim against data, and asks a second question alongside it: once we account for altitude, does the country a coffee comes from still matter?
The dataset and its original source
The data comes from the Coffee Quality Institute (CQI), a nonprofit that trains and licenses the professional coffee graders known as Q Graders. Every row in this file is one coffee sample that a licensed Q Grader evaluated under CQI’s standardized cupping protocol, and the results are published in CQI’s Coffee Quality Program database.
- Original source: Coffee Quality Institute — https://www.coffeeinstitute.org/
- CQI database: https://database.coffeeinstitute.org/
- Where I downloaded the CSV: the TidyTuesday archive for 2020-07-07, https://github.com/rfordatascience/tidytuesday/tree/master/data/2020/2020-07-07, which redistributes a compilation of the CQI database assembled by James LeDoux (https://github.com/jldbc/coffee-quality-database)
TidyTuesday and the coffee-quality-database are repositories that pass the data along. CQI is the body that actually collected and published it.
Variables I use, defined for the reader
The file has 1,339 coffees and 43 variables. These are the ones this project relies on:
Quantitative variables
total_cup_points— the coffee’s overall quality score on a 0–100 scale. In specialty coffee, 80 is the threshold for “specialty grade,” and scores above 85 are rare and commercially valuable.altitude_mean_meters— the midpoint of the altitude range where the coffee was grown, in meters above sea level.aroma,flavor,aftertaste,acidity,body,balance,uniformity,clean_cup,sweetness— the individual attributes the grader scores from 0 to 10, which sum into the total. “Acidity” here is a positive term meaning brightness or liveliness in the cup, not sourness.
Categorical variables
country_of_origin— where the coffee was grown (36 countries).processing_method— how the coffee cherry was pulped and dried: Washed/Wet, Natural/Dry, Semi-washed/Semi-pulped, Pulped natural/honey, or Other.species— Arabica or Robusta.region— a five-level grouping of countries that I construct during cleaning, since the raw file has no continent variable.
What I plan to explore
I want to know whether growing altitude predicts cupping score, and whether that relationship is strong enough to justify how much the industry leans on it. I expect a positive relationship, and I am specifically watching for countries that break the pattern, because those are more interesting than the ones that confirm it.
Setup
# Import with readr::read_csv rather than base read.csv. read_csv returns a tibble,
# parses column types explicitly, and is considerably faster on files this size.
coffee_raw <- readr::read_csv("coffee_ratings.csv")
# Confirm the dimensions match what I expect before touching anything.
dim(coffee_raw)[1] 1339 43Data Cleaning
The raw file needs real work before it can be trusted. I handle each problem in its own step so the reasoning stays visible.
Step 1: Inspect the quality score for impossible values
# A cupping score is bounded between 0 and 100, and anything that reached this
# database should be at least in the 60s. Check the minimum for void records.
range(coffee_raw$total_cup_points)[1] 0.00 90.58# One row scores exactly 0, which is not a possible grade. Look at it directly.
coffee_raw |>
filter(total_cup_points == 0) |>
select(country_of_origin, total_cup_points, aroma, flavor, aftertaste)# A tibble: 1 × 5
country_of_origin total_cup_points aroma flavor aftertaste
<chr> <dbl> <dbl> <dbl> <dbl>
1 Honduras 0 0 0 0That row has zeros across every attribute, so it is a void or failed record rather than a genuinely terrible coffee. It has to go, because leaving it in would drag any mean it touches downward by a meaningless amount.
Step 2: Diagnose the altitude column
# Coffee is not grown above roughly 2,500 m anywhere in the world, so check the
# upper tail of the altitude column for values that cannot be real.
coffee_raw |>
filter(altitude_mean_meters > 3000) |>
select(country_of_origin, altitude, unit_of_measurement, altitude_mean_meters) |>
head(10)# A tibble: 10 × 4
country_of_origin altitude unit_of_measurement altitude_mean_meters
<chr> <chr> <chr> <dbl>
1 Guatemala 3280 m 3280
2 Brazil 11000 metros m 11000
3 Colombia 1800 meters (5900 m 3850
4 Guatemala 3280 m 3280
5 Myanmar 4001 m 4001
6 Guatemala 190164 m 190164
7 Guatemala 3280 m 3280
8 Myanmar 3825 m 3825
9 Nicaragua 1100.00 mosl m 110000
10 Myanmar 3800 m 3800This reveals three separate defects, all created upstream when free-text altitude entries were parsed into a numeric column:
- Lost decimal points. A Guatemalan coffee entered as
1901.64became190164, roughly 22 times the height of Everest. A Nicaraguan entry of1100.00 moslbecame110000. - Feet recorded as meters. Several Myanmar coffees sit at 3,800–4,287 “meters.” Myanmar’s coffee regions top out near 1,500 m, so these are plainly feet that were labelled meters. The repeated Guatemalan value of
3280is the giveaway: 3,280 feet is exactly 1,000 m. - Range text averaged badly. A Colombian entry reading
1800 meters (5900was averaged into 3,850, mixing a metric figure with its own parenthetical conversion to feet.
I cannot recover the true value for each of these individually without guessing, so the defensible move is to mark implausible altitudes as missing rather than invent corrections.
Step 3: Apply the cleaning
coffee <- coffee_raw |>
# Keep only Arabica. The 28 Robusta coffees are graded on a different scale
# and would not be comparable to the Arabica scores.
filter(species == "Arabica") |>
# Drop the void record identified in Step 1 and the single row with no country.
filter(total_cup_points > 0, !is.na(country_of_origin)) |>
# Standardise country names: shorten one official long form, label the two US
# entries readably, and repair a mojibake apostrophe in "Cote d?Ivoire" that
# came from a character-encoding error in the source file.
mutate(country = recode(country_of_origin,
"Tanzania, United Republic Of" = "Tanzania",
"United States (Hawaii)" = "Hawaii (USA)",
"United States (Puerto Rico)" = "Puerto Rico (USA)",
"Cote d?Ivoire" = "Cote d'Ivoire")) |>
# Recode implausible altitudes to NA rather than deleting the rows, so those
# coffees still count in analyses that do not involve altitude.
mutate(altitude_m = if_else(altitude_mean_meters > 3000 | altitude_mean_meters < 100,
NA_real_,
altitude_mean_meters)) |>
# Build the region variable the raw file lacks, so countries can be coloured
# by growing area on the final plot.
mutate(region = case_when(
country %in% c("Ethiopia", "Kenya", "Uganda", "Tanzania", "Malawi",
"Burundi", "Rwanda", "Zambia", "Mauritius",
"Cote d'Ivoire") ~ "Africa",
country %in% c("Colombia", "Brazil", "Peru", "Ecuador") ~ "South America",
country %in% c("Mexico", "Guatemala", "Honduras", "Costa Rica",
"Nicaragua", "El Salvador", "Panama", "Haiti") ~ "Mexico & Central America",
country %in% c("Hawaii (USA)", "Puerto Rico (USA)",
"United States") ~ "United States",
TRUE ~ "Asia & Pacific"))
# Report what survived cleaning and how much altitude data was sacrificed.
cat("Coffees after cleaning:", nrow(coffee), "\n")Coffees after cleaning: 1309 cat("Rows with usable altitude:", sum(!is.na(coffee$altitude_m)), "\n")Rows with usable altitude: 1051 cat("Rows with altitude set to NA:", sum(is.na(coffee$altitude_m)), "\n")Rows with altitude set to NA: 258 Exploratory Analysis
Exploration 1: How are quality scores distributed?
# A histogram is the right first look at a single continuous variable. The dashed
# line marks 80 points, the industry threshold for "specialty grade" coffee.
ggplot(coffee, aes(x = total_cup_points)) +
geom_histogram(binwidth = 0.5, fill = "#6F4E37", color = "white") +
geom_vline(xintercept = 80, linetype = "dashed", color = "#D55E00", linewidth = 0.8) +
annotate("text", x = 79.5, y = 60, label = "Specialty grade threshold",
hjust = 1, size = 3.5, color = "#D55E00") +
labs(title = "Almost every coffee in this database clears specialty grade",
x = "Total Cupping Score (0-100 scale)",
y = "Number of Coffees") +
theme_minimal(base_size = 12)
Scores cluster tightly between 80 and 85. That is a selection effect worth naming early: producers submit coffees to CQI hoping for a good grade, so this is a sample of already-decent coffee, not a sample of all coffee.
Exploration 2: Do all ten attributes actually vary?
# Reshape the ten attribute columns into long format so I can compute, for each
# attribute, the proportion of coffees awarded the maximum score of 10.
ceiling_check <- coffee |>
select(aroma, flavor, aftertaste, acidity, body, balance,
uniformity, clean_cup, sweetness) |>
pivot_longer(everything(), names_to = "attribute", values_to = "score") |>
group_by(attribute) |>
summarize(pct_perfect = 100 * mean(score == 10, na.rm = TRUE), .groups = "drop")
ggplot(ceiling_check, aes(x = pct_perfect, y = reorder(attribute, pct_perfect))) +
geom_col(fill = "#0072B2") +
labs(title = "Three of the ten scored attributes are almost always a perfect 10",
subtitle = "These components cannot distinguish a good coffee from a great one",
x = "Percent of Coffees Scored a Perfect 10",
y = NULL) +
theme_minimal(base_size = 12)
This is the most useful thing exploration turned up. Sweetness, clean cup, and uniformity are scored a flawless 10 for 86–93% of coffees, so they are effectively constants. All the real variation in total_cup_points is carried by aroma, flavor, aftertaste, acidity, body, and balance.
Exploration 3: Does score vary by region?
# Side-by-side boxplots compare a continuous variable across a categorical one,
# showing medians, spread, and outliers together. reorder() sorts regions by
# median score so the comparison reads left to right.
ggplot(coffee, aes(x = reorder(region, total_cup_points, median),
y = total_cup_points)) +
geom_boxplot(fill = "#009E73", alpha = 0.7, outlier.alpha = 0.4) +
coord_flip() +
labs(title = "African coffees score highest as a group",
x = NULL,
y = "Total Cupping Score (0-100 scale)") +
theme_minimal(base_size = 12)
African coffees lead. That points toward altitude, since East African coffee is grown very high, but region and altitude are tangled together and this plot cannot separate them. The final visualization puts both on the same axes.
Final Visualization
# Aggregate to one row per country: median altitude, median score, and the number
# of coffees behind each. Medians resist the influence of a few extreme coffees.
# Countries with fewer than 10 rated coffees are excluded, because a median built
# on three or four samples is too unstable to plot next to Mexico's 231.
country_summary <- coffee |>
filter(!is.na(altitude_m)) |>
group_by(country, region) |>
summarize(n_coffees = n(),
med_alt = median(altitude_m),
med_score = median(total_cup_points),
.groups = "drop") |>
filter(n_coffees >= 10)
cat("Countries plotted:", nrow(country_summary),
"| Coffees represented:", sum(country_summary$n_coffees), "\n")Countries plotted: 17 | Coffees represented: 998 # Correlation at the country level, for reference in the write-up.
cor(country_summary$med_alt, country_summary$med_score)[1] 0.650944# A deliberately chosen palette from the Okabe-Ito colourblind-safe set, one
# colour per growing region. These are intentional, non-default colours.
region_palette <- c("Africa" = "#D55E00",
"South America" = "#0072B2",
"Mexico & Central America" = "#009E73",
"Asia & Pacific" = "#CC79A7",
"United States" = "#E69F00")
final_plot <- ggplot(country_summary, aes(x = med_alt, y = med_score)) +
# Linear trend line first, so it sits behind the points.
geom_smooth(method = "lm", se = TRUE,
color = "grey35", fill = "grey85", linewidth = 0.7) +
# Colour encodes region (categorical); size encodes sample size (quantitative).
geom_point(aes(color = region, size = n_coffees), alpha = 0.9) +
# Repelled labels so every country is identifiable without a lookup table.
geom_text_repel(aes(label = country), size = 3.2, seed = 42,
box.padding = 0.5, max.overlaps = 20,
segment.color = "grey60") +
scale_color_manual(name = "Growing Region", values = region_palette) +
scale_size_continuous(name = "Coffees Rated", range = c(3, 10),
breaks = c(25, 100, 200)) +
scale_x_continuous(breaks = seq(750, 1800, 250)) +
labs(
title = "Higher-Grown Coffee Scores Better, and East Africa Sits at the Top",
subtitle = "Each point is one country, positioned by its median growing altitude and median cupping score",
x = "Median Growing Altitude (meters above sea level)",
y = "Median Cupping Score (0-100 scale)",
caption = "Source: Coffee Quality Institute, Coffee Quality Program database (1,311 Arabica coffees graded by licensed Q Graders).\nCountries shown are those with at least 10 rated coffees and usable altitude data."
) +
# A non-default theme, with further adjustments for readability.
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 15),
plot.subtitle = element_text(color = "grey30", margin = margin(b = 10)),
plot.caption = element_text(color = "grey40", hjust = 0, size = 8.5),
legend.position = "right",
panel.grid.minor = element_blank()
)
final_plot
Essay
How I cleaned the dataset
Cleaning fell into four operations.
Subsetting for comparability. I filtered to the 1,311 Arabica coffees and excluded the 28 Robusta samples. Robusta is graded on a separate CQI scale, so mixing the two would compare scores that do not mean the same thing.
Removing a void record. Examining the range of total_cup_points surfaced one row scoring exactly 0, with zeros across every component attribute. A cupping score of 0 is not attainable, so this is a failed or voided grading rather than a real evaluation, and I dropped it. I also dropped the single row with a missing country_of_origin.
Recoding implausible altitudes as missing. This was the substantive decision. Inspecting the upper tail of altitude_mean_meters exposed three distinct upstream parsing defects: decimal points dropped during import, so 1901.64 became 190164; altitudes reported in feet but labelled meters, which is why several Myanmar coffees claim 3,800–4,287 m and why the value 3280 recurs, since 3,280 feet is exactly 1,000 m; and free-text ranges averaged incorrectly, as with the Colombian entry 1800 meters (5900 that became 3,850. Because I cannot recover any individual true value without guessing, I used if_else() to recode altitudes above 3,000 m or below 100 m to NA instead of imputing a substitute. I set them to missing rather than deleting the rows so those coffees still contribute to analyses that do not involve altitude. This cost me altitude data but preserved honesty about what the file can support.
Standardising labels and deriving a variable. I recoded four country names, including repairing Cote d?Ivoire, where a mojibake artifact had replaced an apostrophe during character encoding. The raw file has no continent or region field, so I derived a five-level region variable with case_when() to group the 36 countries by growing area, which is what lets the final plot use colour to carry a second dimension.
What the visualization represents, and what surprised me
The final plot places each coffee-growing country by two medians: growing altitude on the x-axis and cupping score on the y-axis. Colour distinguishes the growing region, point size shows how many rated coffees stand behind each median, and the grey line is a linear trend with its confidence band. Aggregating to the country level was deliberate — plotting 998 individual coffees produced a cloud too dense to read, while country medians make the structure legible and let each point be labelled.
The main pattern holds up. Higher-altitude countries score better, and the top-right corner belongs to East Africa: Ethiopia at 1,800 m and 85.2 points, Kenya at 1,754 m and 84.5, Uganda at 1,675 m and 83.9. The bottom left is Nicaragua at 1,100 m and 80.9. Ethiopia’s four-point lead over Nicaragua is large in this trade, where a single point can move a lot on price.
The first surprise was how much aggregation flattered the relationship. At the country level the correlation between median altitude and median score is 0.65, which looks impressive. But computed across individual coffees, the same correlation is only 0.24. Aggregating to medians averages away the within-country variation, so the country-level number describes a relationship between country averages, not between coffees. Reading 0.65 as the strength of the altitude effect on any particular coffee would be a mistake — it is close to an ecological correlation, and it made me much less confident in the industry’s “high-grown means good” shorthand than the plot alone suggests. The honest statement is that altitude explains some of the variation between countries and rather little within them.
The second surprise was Tanzania. It sits at 1,600 m, higher than Guatemala or Costa Rica, yet medians only 82.2 points, well below the trend line. Hawaii is a mirror image: Kona coffee grows around 500 m and still medians 82.8, above Mexico at 1,250 m. Two clear counterexamples in seventeen countries is a reminder that altitude is a proxy for the things that actually matter — cooler temperatures, slower cherry maturation, varietal choice, and processing care — and proxies fail when the underlying causes come apart.
The third surprise came from exploration rather than the final plot. Three of the ten scored attributes are nearly useless as discriminators: sweetness, clean cup, and uniformity are awarded a perfect 10 to between 86% and 93% of coffees. A third of the scoring rubric is doing almost no work, and the entire spread in total score is produced by the remaining attributes. Combined with the fact that nearly every coffee here clears the 80-point specialty threshold, this is a database of already-good coffee being separated by fine margins on a handful of dimensions.
What I wanted to show but could not
Hawaii had to leave the altitude analysis. Hawaii contributes 73 coffees, the sixth-largest group in the file, but 71 of them have no altitude recorded. With only two usable values it fell below my ten-coffee threshold and dropped out of the final plot. That is frustrating precisely because Hawaii is the most interesting counterexample I found: low-altitude coffee that scores well. I could describe it in prose but could not show it.
I could not separate altitude from region. These two variables are badly confounded. East African coffee is high-grown, so “Africa scores well” and “high altitude scores well” are nearly the same statement in this dataset, and nothing in the file lets me disentangle them. Doing so properly would need a model with both terms and an interaction, along with varietal and temperature data the file does not include, which is beyond what this project covers.
I wanted to use processing method as a third dimension — my original plan was to colour by processing method rather than region. I abandoned it after checking: median scores across the five methods span barely half a point, from 82.4 for Washed/Wet to 83.0 for Natural/Dry, and 170 coffees have no method recorded at all. Encoding a variable with that little signal would have added colour without adding information, so region earned the colour channel instead.
A time dimension was not available in usable form. harvest_year is free text with entries like 2013/2014, Sept 2009 - April 2010, and mmm, so building a trend over harvest years would have required a parsing effort larger than the rest of the cleaning combined.