tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.cd /app && ls -la && echo "=== OUTPUTS DIR ===" && ls -la outputs/ 2>/dev/null && echo "=== DATA ===" && ls -la data/
total 24 drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 17:46 . drwxr-xr-x 28 root root 4096 Jun 20 18:18 .. -rw-rw---- 1 h2tagent h2tagent 8171 Jun 20 17:44 analysis.R drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:46 data drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:46 outputs === OUTPUTS DIR === total 8 drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:46 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 17:46 .. === DATA === total 598 drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:46 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 17:46 .. -rw-rw---- 1 h2tagent h2tagent 603612 Jun 20 17:44 coffee_ratings.csv
cd /app && cat analysis.R
#!/usr/bin/env Rscript
# Coffee outlier audit -- first draft
# ===================================
# Quick first pass before the deadline. Loaded the CSV, threw mean / sd / IQR
# at each column, used standardised residuals to flag influential lots, and
# took simple country means. Most of this should be good enough but the
# numbers feel off in a few places.
suppressPackageStartupMessages({
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(jsonlite)
library(purrr)
})
DATA_PATH <- "/app/data/coffee_ratings.csv"
OUTPUT_DIR <- "/app/outputs"
dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE)
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
dat <- raw %>% filter(total_cup_points > 0)
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L))
# ---- Altitude
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
q <- quantile(alt_pos, c(0.25, 0.75))
iqr_alt <- q[[2]] - q[[1]]
lo_m <- q[[1]] - 3 * iqr_alt
hi_m <- q[[2]] + 3 * iqr_alt
flag_alt <- !is.na(alt) & (alt < lo_m | alt > hi_m)
n_flag_alt <- sum(flag_alt)
alt_corrected <- rep(NA_real_, nrow(dat))
unit_corrected <- rep(FALSE, nrow(dat))
n_unit_corrected <- 0L
# ---- Defect counts
defect_summary <- function(x) {
q <- quantile(x, c(0.25, 0.75), na.rm = TRUE)
upper <- q[[2]] + 3 * (q[[2]] - q[[1]])
flagged <- !is.na(x) & x > upper
list(upper_threshold = upper, n = sum(flagged), flag_vec = flagged)
}
d_cat1 <- defect_summary(dat$category_one_defects)
d_cat2 <- defect_summary(dat$category_two_defects)
d_quak <- defect_summary(dat$quakers)
flag_defect <- d_cat1$flag_vec | d_cat2$flag_vec | d_quak$flag_vec
# ---- Multivariate grades (every grade column)
g_all <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "uniformity", "clean_cup",
"sweetness", "cupper_points")
G <- as.matrix(dat[, g_all])
G_complete <- complete.cases(G)
mu <- colMeans(G[G_complete, ])
S <- cov(G[G_complete, ])
md2 <- rep(NA_real_, nrow(G))
md2[G_complete] <- mahalanobis(G[G_complete, ], mu, S, tol = 1e-30)
md2_thresh <- qchisq(0.999, df = length(g_all))
flag_mv <- !is.na(md2) & md2 > md2_thresh
n_flag_mv <- sum(flag_mv)
# ---- Influence: standardised residuals
fit_df <- dat[G_complete, ]
fit <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points, data = fit_df)
rstd <- rstandard(fit)
infl_local <- abs(rstd) > 2
n_fit <- length(rstd)
p_fit <- length(coef(fit))
infl_full <- rep(FALSE, nrow(dat))
infl_full[G_complete] <- infl_local
n_high_lev <- sum(infl_full)
n_high_cook <- sum(infl_full)
n_both <- sum(infl_full)
top_idx <- order(-abs(rstd))[1:10]
top10_row_ids <- as.integer(dat$row_id[which(G_complete)[top_idx]])
# ---- Country ranking (raw means)
country_tbl <- dat %>%
filter(!is.na(country_of_origin)) %>%
group_by(country = country_of_origin) %>%
summarise(n = n(),
raw_mean = mean(total_cup_points),
.groups = "drop") %>%
filter(n >= 10) %>%
arrange(desc(raw_mean)) %>%
mutate(raw_rank = rank(-raw_mean, ties.method = "min"),
trimmed_mean = raw_mean,
trimmed_rank = raw_rank,
abs_rank_change = 0L)
n_countries <- nrow(country_tbl)
n_changes_ge_2 <- 0L
max_change <- 0L
# ---- Composite + sensitivity
composite <- flag_alt | flag_mv | infl_full
keep_for_sens <- !composite & G_complete
fit_drop <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points, data = dat[keep_for_sens, ])
coef_full <- unname(coef(fit)["flavor"])
coef_dropped <- unname(coef(fit_drop)["flavor"])
delta_pct <- 100 * (coef_dropped - coef_full) / coef_full
# ---- Flag interactions
A <- flag_alt; M <- flag_mv; D <- flag_defect; I <- infl_full
total_flags <- as.integer(A) + as.integer(M) + as.integer(D) + as.integer(I)
flag_interactions <- list(
n_with_zero_flags = unbox(as.integer(sum(total_flags == 0L))),
…[truncated]cd /app && head -1 data/coffee_ratings.csv | tr ',' '\n' | nl
1 total_cup_points
2 species
3 owner
4 country_of_origin
5 farm_name
6 lot_number
7 mill
8 ico_number
9 company
10 altitude
11 region
12 producer
13 number_of_bags
14 bag_weight
15 in_country_partner
16 harvest_year
17 grading_date
18 owner_1
19 variety
20 processing_method
21 aroma
22 flavor
23 aftertaste
24 acidity
25 body
26 balance
27 uniformity
28 clean_cup
29 sweetness
30 cupper_points
31 moisture
32 category_one_defects
33 quakers
34 color
35 category_two_defects
36 expiration
37 certification_body
38 certification_address
39 certification_contact
40 unit_of_measurement
41 altitude_low_meters
42 altitude_high_meters
43 altitude_mean_meterscd /app && R --version | head -1 && echo "=== Packages ===" && Rscript -e 'for (p in c("dplyr","tidyr","readr","stringr","jsonlite","purrr")) cat(p, as.character(packageVersion(p)), "\n")'R version 4.4.2 (2024-10-31) -- "Pile of Leaves" === Packages === dplyr 1.1.4 tidyr 1.3.1 readr 2.1.5 stringr 1.5.1 jsonlite 1.9.0 purrr 1.0.4
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
cat("n_input =", nrow(raw), "\n")
grades <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
# zero-score lot
allzero <- rowSums(raw[,grades] == 0, na.rm=TRUE) == length(grades)
cat("rows with ALL ten grades == 0:", sum(allzero), "\n")
cat("rows with total_cup_points == 0:", sum(raw$total_cup_points == 0, na.rm=TRUE), "\n")
cat("rows with total_cup_points <= 0:", sum(raw$total_cup_points <= 0, na.rm=TRUE), "\n")
print(raw[allzero, c("total_cup_points","aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")])
'n_input = 1339
rows with ALL ten grades == 0: 1
rows with total_cup_points == 0: 1
rows with total_cup_points <= 0: 1
# A tibble: 1 × 11
total_cup_points aroma flavor aftertaste acidity body balance uniformity
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0 0 0 0 0 0 0
# ℹ 3 more variables: clean_cup <dbl>, sweetness <dbl>, cupper_points <dbl>cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
grades <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
dat <- raw %>% filter(total_cup_points > 0)
cat("n_after_drop =", nrow(dat), "\n\n")
cat("=== Grade column summaries (note near-constants) ===\n")
for (g in grades) {
x <- dat[[g]]
cat(sprintf("%-14s mean=%.3f sd=%.4f min=%.2f max=%.2f\n", g, mean(x,na.rm=T), sd(x,na.rm=T), min(x,na.rm=T), max(x,na.rm=T)))
}
cat("\n=== uniformity/clean_cup/sweetness value tables (why near-constant) ===\n")
for (g in c("uniformity","clean_cup","sweetness")) { cat(g, ":\n"); print(table(dat[[g]])) }
'n_after_drop = 1338 === Grade column summaries (note near-constants) === aroma mean=7.572 sd=0.3159 min=5.08 max=8.75 flavor mean=7.526 sd=0.3414 min=6.08 max=8.83 aftertaste mean=7.407 sd=0.3503 min=6.17 max=8.67 acidity mean=7.541 sd=0.3192 min=5.25 max=8.75 body mean=7.523 sd=0.3078 min=5.08 max=8.58 balance mean=7.524 sd=0.3536 min=5.25 max=8.75 uniformity mean=9.842 sd=0.4852 min=6.00 max=10.00 clean_cup mean=9.842 sd=0.7153 min=0.00 max=10.00 sweetness mean=9.864 sd=0.5542 min=1.33 max=10.00 cupper_points mean=7.509 sd=0.4268 min=5.17 max=10.00 === uniformity/clean_cup/sweetness value tables (why near-constant) === uniformity : 6 6.67 7.33 8 8.67 9 9.33 9.5 10 3 7 2 25 31 1 116 1 1152 clean_cup : 0 1.33 2.67 5.33 6 6.67 7.33 8 8.67 9.33 10 1 1 2 3 6 13 3 13 16 61 1219 sweetness : 1.33 6 6.67 7.08 7.42 7.5 7.58 7.67 7.75 7.83 7.92 8 8.42 8.67 9.33 10 1 3 8 1 3 1 5 2 7 2 2 11 1 12 61 1218
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
for (g in c("category_one_defects","category_two_defects","quakers")) {
x <- dat[[g]]
cat("===", g, "=== n_NA =", sum(is.na(x)), "\n")
cat(sprintf(" pct zero = %.1f%%\n", 100*mean(x==0,na.rm=T)))
print(quantile(x, c(.5,.75,.9,.95,.98,.99,1), na.rm=TRUE))
q <- quantile(x, c(.25,.75), na.rm=TRUE); cat(" raw Tukey k=3 upper =", q[[2]]+3*(q[[2]]-q[[1]]), " (Q1,Q3=",q[[1]],q[[2]],")\n\n")
}
'=== category_one_defects === n_NA = 0 pct zero = 84.9% 50% 75% 90% 95% 98% 99% 100% 0.00 0.00 1.00 2.00 5.00 7.63 63.00 raw Tukey k=3 upper = 0 (Q1,Q3= 0 0 ) === category_two_defects === n_NA = 0 pct zero = 27.9% 50% 75% 90% 95% 98% 99% 100% 2.00 4.00 9.00 13.00 20.00 26.63 55.00 raw Tukey k=3 upper = 16 (Q1,Q3= 0 4 ) === quakers === n_NA = 1 pct zero = 93.0% 50% 75% 90% 95% 98% 99% 100% 0 0 0 1 2 5 11 raw Tukey k=3 upper = 0 (Q1,Q3= 0 0 )
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
n <- nrow(dat)
for (g in c("category_one_defects","category_two_defects","quakers")) {
x <- dat[[g]]; xnn <- x[!is.na(x)]
pos <- xnn[xnn>0]
cat("===", g, "=== (n nonzero =", length(pos), ")\n")
# Option A: Tukey k=1.5 & k=3 on positive subset
for (k in c(1.5,3)) {
q <- quantile(pos, c(.25,.75)); up <- q[[2]]+k*(q[[2]]-q[[1]])
cat(sprintf(" Tukey-on-positives k=%.1f: upper=%.2f n_flag(>up)=%d (%.2f%%)\n", k, up, sum(xnn>up), 100*sum(xnn>up)/n))
}
# Option C: quantiles of full column
for (p in c(.98,.99)) {
th <- quantile(xnn, p); cat(sprintf(" q%.0f full: thr=%.2f n_flag(>thr)=%d (%.2f%%)\n", 100*p, th, sum(xnn>th), 100*sum(xnn>th)/n))
}
cat("\n")
}
'=== category_one_defects === (n nonzero = 202 ) Tukey-on-positives k=1.5: upper=6.00 n_flag(>up)=17 (1.27%) Tukey-on-positives k=3.0: upper=9.00 n_flag(>up)=12 (0.90%) q98 full: thr=5.00 n_flag(>thr)=20 (1.49%) q99 full: thr=7.63 n_flag(>thr)=14 (1.05%) === category_two_defects === (n nonzero = 965 ) Tukey-on-positives k=1.5: upper=12.00 n_flag(>up)=73 (5.46%) Tukey-on-positives k=3.0: upper=18.00 n_flag(>up)=34 (2.54%) q98 full: thr=20.00 n_flag(>thr)=23 (1.72%) q99 full: thr=26.63 n_flag(>thr)=14 (1.05%) === quakers === (n nonzero = 94 ) Tukey-on-positives k=1.5: upper=6.00 n_flag(>up)=6 (0.45%) Tukey-on-positives k=3.0: upper=9.00 n_flag(>up)=1 (0.07%) q98 full: thr=2.00 n_flag(>thr)=25 (1.87%) q99 full: thr=5.00 n_flag(>thr)=10 (0.75%)
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
alt <- dat$altitude_mean_meters
cat("altitude_mean_meters: n_NA =", sum(is.na(alt)), " n_pos =", sum(alt>0,na.rm=T), " n_zero =", sum(alt==0,na.rm=T),"\n")
cat("skewness check: mean=",mean(alt,na.rm=T)," median=",median(alt,na.rm=T),"\n")
print(quantile(alt, c(0,.25,.5,.75,.9,.99,1), na.rm=TRUE))
cat("\n=== Raw-meter Tukey k=3 (the DRAFT approach) ===\n")
ap <- alt[!is.na(alt)&alt>0]; q<-quantile(ap,c(.25,.75)); iqr<-q[[2]]-q[[1]]
cat("raw lo=",q[[1]]-3*iqr," hi=",q[[2]]+3*iqr,"\n")
cat("\n=== log10-scale Tukey k=3 (REQUIRED approach) ===\n")
L <- log10(ap); qL<-quantile(L,c(.25,.75)); iqrL<-qL[[2]]-qL[[1]]
loL<-qL[[1]]-3*iqrL; hiL<-qL[[2]]+3*iqrL
cat("log10 Q1,Q3 =",qL[[1]],qL[[3]]," IQR=",iqrL,"\n")
cat("bounds log10: lo=",loL," hi=",hiL,"\n")
cat("bounds meters: lo=",10^loL," hi=",10^hiL,"\n")
flag <- !is.na(alt) & alt>0 & (log10(alt) < loL | log10(alt) > hiL)
flagna <- !is.na(alt) & alt<=0
cat("n flagged (outside fence, pos only) =", sum(flag), "\n")
cat("rows with alt<=0:", sum(alt<=0,na.rm=T), "\n")
cat("\nTop extreme altitudes:\n")
print(sort(unique(ap[ap>10^hiL])))
print(sort(unique(ap[ap<10^loL])))
'altitude_mean_meters: n_NA = 230 n_pos = 1108 n_zero = 0
skewness check: mean= 1775.369 median= 1310.64
0% 25% 50% 75% 90% 99% 100%
1.00 1100.00 1310.64 1600.00 1799.30 3280.00 190164.00
=== Raw-meter Tukey k=3 (the DRAFT approach) ===
raw lo= -400 hi= 3100
=== log10-scale Tukey k=3 (REQUIRED approach) ===
[stderr]
Error in qL[[3]] : subscript out of bounds
Calls: cat
Execution haltedcd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
alt <- dat$altitude_mean_meters
ap <- alt[!is.na(alt)&alt>0]
L <- log10(ap); qL<-quantile(L,c(.25,.75)); iqrL<-qL[[2]]-qL[[1]]
loL<-qL[[1]]-3*iqrL; hiL<-qL[[2]]+3*iqrL
cat(sprintf("log10 Q1=%.4f Q3=%.4f IQR=%.4f\n", qL[[1]],qL[[2]],iqrL))
cat(sprintf("bounds meters: lo=%.4f hi=%.4f\n", 10^loL, 10^hiL))
flag <- !is.na(alt) & alt>0 & (log10(alt) < loL | log10(alt) > hiL)
cat("n flagged =", sum(flag), "\n\n")
df <- dat[flag, c("altitude","altitude_mean_meters","unit_of_measurement")]
print(df, n=40)
'log10 Q1=3.0414 Q3=3.2041 IQR=0.1627 bounds meters: lo=357.4463 hi=4923.8167 n flagged = 51 # A tibble: 51 × 3 altitude altitude_mean_meters unit_of_measurement <chr> <dbl> <chr> 1 1.2 12 m 2 1.2 12 m 3 -1 1 m 4 350 350 m 5 170 170 m 6 approx. 1'500 metres a.s.l. 250. m 7 1 1 m 8 175 175 m 9 16.8 168 m 10 approx.1'500 metres a.s.l. 250. m 11 1 1 m 12 350 350 m 13 11000 metros 11000 m 14 518 158. ft 15 1.65 165 m 16 1 1 m 17 200 m 200 m 18 518 158. ft 19 1 1 m 20 200 m 200 m 21 1 1 m 22 1 1 m 23 1 1 m 24 1.3 13 m 25 518 158. ft 26 1 1 m 27 1 1 m 28 1 1 m 29 150 m 150 m 30 190164 190164 m 31 1.3 13 m 32 160-200m 180 m 33 1.25 125 m 34 110公尺 110 m 35 1 1 m 36 200 200 m 37 300 m 300 m 38 100 100 m 39 1100.00 mosl 110000 m 40 1.85 185 m # ℹ 11 more rows
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(stringr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
alt <- dat$altitude_mean_meters
ap <- alt[!is.na(alt)&alt>0]
L <- log10(ap); qL<-quantile(L,c(.25,.75)); iqrL<-qL[[2]]-qL[[1]]
lo<-10^(qL[[1]]-3*iqrL); hi<-10^(qL[[2]]+3*iqrL)
flag <- !is.na(alt) & alt>0 & (log10(alt) < log10(lo) | log10(alt) > log10(hi))
first_tok <- function(s){ m <- str_extract(s, "[0-9]+\\.?[0-9]*"); as.numeric(m) }
corr <- rep(NA_real_, nrow(dat))
for (i in which(flag)) {
tok <- first_tok(dat$altitude[i])
if (is.na(tok)) next
for (cand in c(tok/10, tok/100, tok)) {
if (!is.na(cand) && cand>=lo && cand<=hi) { corr[i] <- cand; break }
}
}
cat(sprintf("fence lo=%.4f hi=%.4f n_flag=%d n_corrected=%d\n", lo, hi, sum(flag), sum(!is.na(corr))))
res <- dat[flag, c("altitude","altitude_mean_meters")]
res$first_token <- sapply(dat$altitude[flag], first_tok)
res$corrected <- corr[flag]
print(res, n=60)
'fence lo=357.4463 hi=4923.8167 n_flag=51 n_corrected=7 # A tibble: 51 × 4 altitude altitude_mean_meters first_token corrected <chr> <dbl> <dbl> <dbl> 1 1.2 12 1.2 NA 2 1.2 12 1.2 NA 3 -1 1 1 NA 4 350 350 350 NA 5 170 170 170 NA 6 approx. 1'500 metres a.s.l. 250. 1 NA 7 1 1 1 NA 8 175 175 175 NA 9 16.8 168 16.8 NA 10 approx.1'500 metres a.s.l. 250. 1 NA 11 1 1 1 NA 12 350 350 350 NA 13 11000 metros 11000 11000 1100 14 518 158. 518 518 15 1.65 165 1.65 NA 16 1 1 1 NA 17 200 m 200 200 NA 18 518 158. 518 518 19 1 1 1 NA 20 200 m 200 200 NA 21 1 1 1 NA 22 1 1 1 NA 23 1 1 1 NA 24 1.3 13 1.3 NA 25 518 158. 518 518 26 1 1 1 NA 27 1 1 1 NA 28 1 1 1 NA 29 150 m 150 150 NA 30 190164 190164 190164 1902. 31 1.3 13 1.3 NA 32 160-200m 180 160 NA 33 1.25 125 1.25 NA 34 110公尺 110 110 NA 35 1 1 1 NA 36 200 200 200 NA 37 300 m 300 300 NA 38 100 100 100 NA 39 1100.00 mosl 110000 1100 1100 40 1.85 185 1.85 NA 41 50 50 50 NA 42 280 280 280 NA 43 200 m 200 200 NA 44 150 m 150 150 NA 45 250 m 250 250 NA 46 250 m 250 250 NA 47 1901.64 190164 1902. 1902. 48 200 200 200 NA 49 12oo 12 12 NA 50 ~350m 350 350 NA 51 40 40 40 NA
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G <- as.matrix(dat[,g7]); cc <- complete.cases(G)
cat("complete cases on 7 grades:", sum(cc), "of", nrow(dat), "\n")
mu<-colMeans(G[cc,]); S<-cov(G[cc,])
md2<-rep(NA,nrow(dat)); md2[cc]<-mahalanobis(G[cc,],mu,S)
thr<-qchisq(0.999,df=7)
cat("threshold chisq(.999,df=7)=",thr," n_flagged=",sum(!is.na(md2)&md2>thr),"\n\n")
# country
ct <- dat %>% filter(!is.na(country_of_origin)) %>% group_by(country_of_origin) %>% summarise(n=n()) %>% filter(n>=10) %>% arrange(desc(n))
cat("countries with >=10 lots:", nrow(ct), "\n"); print(ct, n=40)
cat("\nrows with NA country:", sum(is.na(dat$country_of_origin)), "\n")
'complete cases on 7 grades: 1338 of 1338 threshold chisq(.999,df=7)= 24.32189 n_flagged= 40 countries with >=10 lots: 21 # A tibble: 21 × 2 country_of_origin n <chr> <int> 1 Mexico 236 2 Colombia 183 3 Guatemala 181 4 Brazil 132 5 Taiwan 75 6 United States (Hawaii) 73 7 Honduras 52 8 Costa Rica 51 9 Ethiopia 44 10 Tanzania, United Republic Of 40 11 Uganda 36 12 Thailand 32 13 Nicaragua 26 14 Kenya 25 15 El Salvador 21 16 Indonesia 20 17 China 16 18 India 14 19 Malawi 11 20 Peru 10 21 United States 10 rows with NA country: 1
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
g10 <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
s <- rowSums(dat[,g10])
cat("max |total - sum(10 grades)| =", max(abs(dat$total_cup_points - s)), "\n")
# quick AIC of the 3 models
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
cc <- complete.cases(dat[,g7]); fr <- dat[cc,]
m1 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points, data=fr)
# altitude
altv <- fr$altitude_mean_meters; le <- log10(altv); le[!is.finite(le)] <- mean(le[is.finite(le)]); fr$log_elev <- le
m2 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+log_elev, data=fr)
for (d in c("category_one_defects","category_two_defects","quakers")){x<-fr[[d]];x[is.na(x)]<-mean(x,na.rm=T);fr[[paste0("l1p_",d)]]<-log1p(x)}
m3 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1p_category_one_defects+l1p_category_two_defects+l1p_quakers, data=fr)
cat(sprintf("AIC grades_only=%.3f grades_altitude=%.3f grades_defects=%.3f\n", AIC(m1),AIC(m2),AIC(m3)))
cat(sprintf("n=%d resid df m1=%d\n", nrow(fr), m1$df.residual))
'max |total - sum(10 grades)| = 0.5 AIC grades_only=4459.362 grades_altitude=4461.275 grades_defects=4418.680 n=1338 resid df m1=1330
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
for (tr in c(0.1, 0.2)) {
ct <- dat %>% filter(!is.na(country_of_origin)) %>% group_by(country=country_of_origin) %>%
summarise(n=n(), raw_mean=mean(total_cup_points), trimmed_mean=mean(total_cup_points,trim=tr), .groups="drop") %>% filter(n>=10)
ct <- ct %>% mutate(raw_rank=rank(-raw_mean,ties.method="min"), trimmed_rank=rank(-trimmed_mean,ties.method="min"), chg=abs(raw_rank-trimmed_rank))
cat(sprintf("trim=%.1f: n_changes>=2 = %d, max_change = %d\n", tr, sum(ct$chg>=2), max(ct$chg)))
}
'trim=0.1: n_changes>=2 = 2, max_change = 2 trim=0.2: n_changes>=2 = 2, max_change = 2
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
fr <- dat # all complete on g7
for (d in c("category_one_defects","category_two_defects","quakers")){x<-fr[[d]];x[is.na(x)]<-mean(x,na.rm=T);fr[[paste0("l1p_",d)]]<-log1p(x)}
m3 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1p_category_one_defects+l1p_category_two_defects+l1p_quakers, data=fr)
n <- nobs(m3); p <- length(coef(m3))
lev <- hatvalues(m3); ck <- cooks.distance(m3)
lev_t <- 2*p/n; ck_t <- 4/n
cat(sprintf("n=%d p=%d lev_thr=%.6f cook_thr=%.6f\n", n,p,lev_t,ck_t))
cat(sprintf("n_high_lev=%d n_high_cook=%d n_both=%d\n", sum(lev>lev_t), sum(ck>ck_t), sum(lev>lev_t & ck>ck_t)))
top10 <- order(-ck)[1:10]
cat("top10 cook row_ids (0-indexed):", paste(top10-1, collapse=","), "\n")
cat("top10 cook values:", paste(round(ck[top10],4),collapse=","), "\n")
'n=1338 p=11 lev_thr=0.016442 cook_thr=0.002990 n_high_lev=110 n_high_cook=66 n_both=35 top10 cook row_ids (0-indexed): 1309,1308,1333,1295,1306,1260,1307,1302,1291,1259 top10 cook values: 0.3795,0.15,0.0979,0.0936,0.0765,0.0531,0.0472,0.0457,0.0361,0.0313
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(stringr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
# altitude flag
alt<-dat$altitude_mean_meters; ap<-alt[!is.na(alt)&alt>0]; L<-log10(ap);qL<-quantile(L,c(.25,.75));iqrL<-qL[[2]]-qL[[1]]
lo<-10^(qL[[1]]-3*iqrL);hi<-10^(qL[[2]]+3*iqrL)
flag_alt<-!is.na(alt)&alt>0&(alt<lo|alt>hi)
# mv flag
G<-as.matrix(dat[,g7]);cc<-complete.cases(G);mu<-colMeans(G[cc,]);S<-cov(G[cc,]);md2<-rep(NA,nrow(dat));md2[cc]<-mahalanobis(G[cc,],mu,S);flag_mv<-!is.na(md2)&md2>qchisq(.999,7)
# influence (grades_defects)
fr<-dat;for(d in c("category_one_defects","category_two_defects","quakers")){x<-fr[[d]];x[is.na(x)]<-mean(x,na.rm=T);fr[[paste0("l1p_",d)]]<-log1p(x)}
m3<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1p_category_one_defects+l1p_category_two_defects+l1p_quakers,data=fr)
n<-nobs(m3);p<-length(coef(m3));lev<-hatvalues(m3)>2*p/n;ck<-cooks.distance(m3)>4/n;infl<-lev&ck
composite<-flag_alt|flag_mv|infl
cat("n composite =",sum(composite)," (alt",sum(flag_alt),"mv",sum(flag_mv),"infl",sum(infl),")\n")
m_full<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat)
m_drop<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat[!composite,])
cf<-coef(m_full)["flavor"];cd<-coef(m_drop)["flavor"]
cat(sprintf("coef_full=%.6f coef_drop=%.6f delta_pct=%.4f\n",cf,cd,100*(cd-cf)/cf))
'n composite = 105 (alt 51 mv 40 infl 35 ) coef_full=2.085398 coef_drop=1.652728 delta_pct=-20.7476
command: create
path: /app/analysis.R
file_text: #!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-matched rewrite
# =====================================================
# Each column is audited with a method that matches its distributional shape,
# replacing the draft's off-the-shelf summaries. See audit_memo.md for the
# rationale behind every rejection. This script alone reproduces all artifacts
# from a clean /app/outputs/ directory.
suppressPackageStartupMessages({
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(jsonlite)
library(purrr)
})
DATA_PATH <- "/app/data/coffee_ratings.csv"
OUTPUT_DIR <- "/app/outputs"
dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE)
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
# ---------------------------------------------------------------------------
# 1. Sentinel drop -- one withdrawn lot has every grade recorded as 0.
# ---------------------------------------------------------------------------
GRADES10 <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"uniformity", "clean_cup", "sweetness", "cupper_points")
all_zero <- rowSums(raw[, GRADES10] == 0, na.rm = TRUE) == length(GRADES10)
dat <- raw[!all_zero, , drop = FALSE]
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L))
N <- nrow(dat)
# ---------------------------------------------------------------------------
# 2. Altitude -- log10-scale Tukey fence (k = 3); meter-scale fences mislead
# on this right-skewed column. Recover decimal-displacement unit slips.
# ---------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
Lq <- quantile(log10(alt_pos), c(0.25, 0.75))
L_iqr <- Lq[[2]] - Lq[[1]]
alt_lo_m <- 10^(Lq[[1]] - 3 * L_iqr)
alt_hi_m <- 10^(Lq[[2]] + 3 * L_iqr)
flag_alt <- !is.na(alt) & alt > 0 & (alt < alt_lo_m | alt > alt_hi_m)
# First numeric token of the raw altitude string (digits + optional decimals;
# thousands-separators/units are intentionally not treated as part of it).
first_token <- function(s) as.numeric(str_extract(s, "[0-9]+\\.?[0-9]*"))
alt_corrected <- rep(NA_real_, N)
for (i in which(flag_alt)) {
tok <- first_token(dat$altitude[i])
if (is.na(tok)) next
for (cand in c(tok / 10, tok / 100, tok)) { # plausible power-of-ten fixes
if (!is.na(cand) && cand >= alt_lo_m && cand <= alt_hi_m) {
alt_corrected[i] <- cand
break
}
}
}
n_flag_alt <- sum(flag_alt)
n_unit_corrected <- sum(!is.na(alt_corrected))
# ---------------------------------------------------------------------------
# 3. Defects -- mass-at-zero counts with a sparse upper tail. Raw Tukey fences
# collapse to 0 (Q1 = Q3 = 0) and flag every nonzero row, so use a per-column
# empirical 98th-percentile upper threshold; flag counts strictly above it.
# ---------------------------------------------------------------------------
DEFECT_COLS <- c("category_one_defects", "category_two_defects", "quakers")
defect_flag_one <- function(x) {
thr <- unname(quantile(x, 0.98, na.rm = TRUE, type = 7))
list(upper_threshold = thr,
flag_vec = !is.na(x) & x > thr,
n = sum(!is.na(x) & x > thr))
}
defect_res <- lapply(DEFECT_COLS, function(cn) defect_flag_one(dat[[cn]]))
names(defect_res) <- DEFECT_COLS
flag_defect <- Reduce(`|`, lapply(defect_res, `[[`, "flag_vec"))
# ---------------------------------------------------------------------------
# 4. Multivariate grades -- classical Mahalanobis on the SEVEN primary SCA
# attributes only (the cleanliness columns are near-constant and inflate /
# destabilise the covariance). Chi-square upper tail at 99.9%, df = 7.
# ---------------------------------------------------------------------------
MV_COLS <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"cupper_points")
G <-
…[truncated]Existing file was successfully overwritten at /app/analysis.R
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R
Audit complete: n_input=1339, n_after=1338, selected=grades_defects, composite=105
cd /app && ls -la outputs/ && echo "=== JSON ===" && cat outputs/outlier_report.json
total 226
drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:22 .
drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:22 ..
-rw-r--r-- 1 h2tagent h2tagent 137312 Jun 20 18:22 coffee_ols_selected.rds
-rw-r--r-- 1 h2tagent h2tagent 78269 Jun 20 18:22 outlier_flags.csv
-rw-r--r-- 1 h2tagent h2tagent 6550 Jun 20 18:22 outlier_report.json
=== JSON ===
{
"n_input": 1339,
"n_after_zero_score_drop": 1338,
"methods": {
"altitude": "log10-scale Tukey/IQR fence (k=3) on positive altitude_mean_meters, bounds back-transformed to meters; flagged rows' first numeric token tested for power-of-ten recovery (/10, /100, as-is) into the fence",
"defects": "Per-column empirical 98th-percentile upper threshold on the mass-at-zero counts (raw IQR fences collapse since Q1=Q3=0); flag counts strictly above threshold",
"multivariate_grades": "Classical Mahalanobis distance (sample mean/covariance) on the seven primary SCA attributes only vs chi-square(0.999, df=7) upper tail",
"influence": "AIC-selected OLS of total_cup_points; flag high leverage (>2p/n) AND high Cook's D (>4/n), n=fit sample size, p=coefficients incl. intercept",
"country_ranking": "Raw mean vs 20% trimmed (robust) mean of total_cup_points for countries with >=10 lots; rank both descending and compare positions"
},
"altitude": {
"lower_bound_m": 357.4462890625,
"upper_bound_m": 4923.8166791886,
"n_flagged": 51,
"n_unit_corrected": 7
},
"defects": {
"category_one_defects": {
"upper_threshold": 5,
"n_flagged": 20
},
"category_two_defects": {
"upper_threshold": 20,
"n_flagged": 23
},
"quakers": {
"upper_threshold": 2,
"n_flagged": 25
}
},
"multivariate_grades": {
"columns_used": ["aroma", "flavor", "aftertaste", "acidity", "body", "balance", "cupper_points"],
"threshold": 24.3218863479,
"n_flagged": 40
},
"influence": {
"selected_model": "grades_defects",
"aic_grades_only": 4459.3615159695,
"aic_grades_altitude": 4461.292329589,
"aic_grades_defects": 4418.6800742954,
"leverage_threshold": 0.01644245142,
"cook_threshold": 0.0029895366218,
"n_high_leverage": 110,
"n_high_cook": 66,
"n_both": 35,
"top10_cook_row_ids": [1309, 1308, 1333, 1295, 1306, 1260, 1307, 1302, 1291, 1259]
},
"country_ranking": {
"n_countries_evaluated": 21,
"n_rank_changes_ge_2": 2,
"max_abs_rank_change": 2,
"per_country": [
{
"country": "Ethiopia",
"n": 44,
"raw_mean": 85.484091,
"trimmed_mean": 85.397857,
"raw_rank": 1,
"trimmed_rank": 1
},
{
"country": "United States",
"n": 10,
"raw_mean": 84.433,
"trimmed_mean": 85.331667,
"raw_rank": 2,
"trimmed_rank": 2
},
{
"country": "Kenya",
"n": 25,
"raw_mean": 84.3096,
"trimmed_mean": 84.544,
"raw_rank": 3,
"trimmed_rank": 3
},
{
"country": "Uganda",
"n": 36,
"raw_mean": 83.451944,
"trimmed_mean": 83.447727,
"raw_rank": 4,
"trimmed_rank": 4
},
{
"country": "Colombia",
"n": 183,
"raw_mean": 83.106557,
"trimmed_mean": 83.267477,
"raw_rank": 5,
"trimmed_rank": 5
},
{
"country": "El Salvador",
"n": 21,
"raw_mean": 83.052857,
"trimmed_mean": 83.059231,
"raw_rank": 6,
"trimmed_rank": 7
},
{
"country": "China",
"n": 16,
"raw_mean": 82.9275,
"trimmed_mean": 83.267,
"raw_rank": 7,
"trimmed_rank": 6
},
{
"country": "Costa Rica",
"n": 51,
"raw_mean": 82.78902,
"trimmed_mean": 83.058065,
"raw_rank": 8,
"trimmed_rank": 8
},
{
"country": "Thailand",
"n": 32,
"raw_mean": 82.57375,
"trimmed_mean": 8
…[truncated]cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr)})
fit <- readRDS("outputs/coffee_ols_selected.rds")
csv <- read_csv("outputs/outlier_flags.csv", show_col_types = FALSE)
cat("lm class:", class(fit), " nobs:", nobs(fit), " coefs:", length(coef(fit)), "\n")
cat("formula:", deparse(formula(fit)), "\n\n")
n <- nobs(fit); p <- length(coef(fit))
lev <- hatvalues(fit); ck <- cooks.distance(fit)
hl <- lev > 2*p/n; hc <- ck > 4/n
# fit rows align with all post-drop rows (reg_idx = all) -> compare directly
cat("CSV rows:", nrow(csv), " fit rows:", length(lev), "\n")
cat("high_leverage match CSV:", all(hl == csv$high_leverage_flag), "\n")
cat("high_cook match CSV:", all(hc == csv$high_cook_flag), "\n")
cat("n_high_lev:", sum(hl), " n_high_cook:", sum(hc), " n_both:", sum(hl&hc), "\n")
# composite reconciliation
comp <- csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag | (csv$high_leverage_flag & csv$high_cook_flag)
cat("composite matches CSV column:", all(comp == csv$composite_drop_flag), "\n")
# top10 cook
cat("top10 cook row_ids:", paste(csv$row_id[order(-ck)[1:10]], collapse=","), "\n")
'lm class: lm nobs: 1338 coefs: 11 formula: total_cup_points ~ aroma + flavor + aftertaste + acidity + body + balance + cupper_points + l1p_category_one_defects + l1p_category_two_defects + l1p_quakers CSV rows: 1338 fit rows: 1338 high_leverage match CSV: TRUE high_cook match CSV: TRUE n_high_lev: 110 n_high_cook: 66 n_both: 35 composite matches CSV column: TRUE top10 cook row_ids: 1309,1308,1333,1295,1306,1260,1307,1302,1291,1259
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types = FALSE)
j <- fromJSON("outputs/outlier_report.json")
cat("=== CSV columns ===\n"); print(names(csv))
cat("\n=== altitude_corrected_m non-NA rows ===\n")
print(csv[!is.na(csv$altitude_corrected_m), c("row_id","altitude_corrected_m","altitude_outlier_flag")])
cat("all corrected rows are altitude-flagged:", all(csv$altitude_outlier_flag[!is.na(csv$altitude_corrected_m)]), "\n\n")
cat("=== JSON vs CSV reconciliation ===\n")
cat("alt flagged:", j$altitude$n_flagged, "==", sum(csv$altitude_outlier_flag), "\n")
cat("alt corrected:", j$altitude$n_unit_corrected, "==", sum(!is.na(csv$altitude_corrected_m)), "\n")
cat("defect flagged (any):", sum(csv$defect_outlier_flag), "\n")
cat("mv flagged:", j$multivariate_grades$n_flagged, "==", sum(csv$multivariate_grade_outlier_flag), "\n")
cat("n_both infl:", j$influence$n_both, "==", sum(csv$high_leverage_flag & csv$high_cook_flag), "\n")
cat("composite:", sum(csv$composite_drop_flag), "\n")
# interactions sum check
fi <- j$flag_interactions
cat("interaction buckets sum:", fi$n_with_zero_flags+fi$n_with_one_flag+fi$n_with_two_flags+fi$n_with_three_flags+fi$n_with_four_flags, "== 1338\n")
cat("types:", sapply(csv[,c("altitude_outlier_flag","altitude_corrected_m","defect_outlier_flag")], class), "\n")
'=== CSV columns === [1] "row_id" "country_of_origin" [3] "total_cup_points" "altitude_outlier_flag" [5] "altitude_corrected_m" "defect_outlier_flag" [7] "multivariate_grade_outlier_flag" "high_leverage_flag" [9] "high_cook_flag" "composite_drop_flag" === altitude_corrected_m non-NA rows === # A tibble: 7 × 3 row_id altitude_corrected_m altitude_outlier_flag <dbl> <dbl> <lgl> 1 543 1100 TRUE 2 584 518 TRUE 3 677 518 TRUE 4 833 518 TRUE 5 896 1902. TRUE 6 1040 1100 TRUE 7 1144 1902. TRUE all corrected rows are altitude-flagged: TRUE === JSON vs CSV reconciliation === alt flagged: 51 == 51 alt corrected: 7 == 7 defect flagged (any): 63 mv flagged: 40 == 40 n_both infl: 35 == 35 composite: 105 interaction buckets sum: 1338 == 1338 types: logical numeric logical
cd /app && cp outputs/outlier_report.json /tmp/run1.json && cp outputs/outlier_flags.csv /tmp/run1.csv && rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1 && echo "JSON identical:" && diff -q /tmp/run1.json outputs/outlier_report.json && echo "CSV identical:" && diff -q /tmp/run1.csv outputs/outlier_flags.csv && echo "Both reproduced byte-for-byte."
JSON identical: CSV identical: Both reproduced byte-for-byte.
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
# Draft raw-meter altitude fence
ap<-dat$altitude_mean_meters; ap<-ap[!is.na(ap)&ap>0]; q<-quantile(ap,c(.25,.75))
cat("DRAFT raw-meter fence: [",q[[1]]-3*(q[[2]]-q[[1]]),",",q[[2]]+3*(q[[2]]-q[[1]]),"] -> flags",
sum(dat$altitude_mean_meters < (q[[2]]+3*(q[[2]]-q[[1]])) & FALSE, na.rm=T), "\n")
draftflag <- !is.na(dat$altitude_mean_meters) & (dat$altitude_mean_meters < (q[[1]]-3*(q[[2]]-q[[1]])) | dat$altitude_mean_meters > (q[[2]]+3*(q[[2]]-q[[1]])))
cat(" draft raw fence flags n =", sum(draftflag), "(misses sub-358m slips, keeps lower bound negative)\n")
# Draft raw-defect Tukey flags
for (g in c("category_one_defects","category_two_defects","quakers")){
x<-dat[[g]]; q<-quantile(x,c(.25,.75),na.rm=T); up<-q[[2]]+3*(q[[2]]-q[[1]])
cat(sprintf(" DRAFT raw Tukey %s: upper=%.1f flags=%d (%.1f%%)\n",g,up,sum(!is.na(x)&x>up),100*sum(!is.na(x)&x>up)/nrow(dat)))
}
# Draft Mahalanobis on 10 cols incl near-constants
g10<-c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
G<-as.matrix(dat[,g10]);cc<-complete.cases(G);mu<-colMeans(G[cc,]);S<-cov(G[cc,])
md2<-mahalanobis(G[cc,],mu,S,tol=1e-30);cat(" DRAFT MV(10 cols) flags n=",sum(md2>qchisq(.999,10)),"(cond # of cov:",round(kappa(S)),")\n")
cat(" NEW MV(7 cols) cond # of cov:",round(kappa(cov(as.matrix(dat[,g10[c(1:6,10)]])[cc,]))),"\n")
# abs(rstandard)>2 count (draft influence)
fit<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat)
cat(" DRAFT abs(rstandard)>2 flags n=",sum(abs(rstandard(fit))>2),"(~",round(100*sum(abs(rstandard(fit))>2)/nrow(dat),1),"% -- expected ~5% by construction)\n")
'DRAFT raw-meter fence: [ -400 , 3100 ] -> flags 0 draft raw fence flags n = 16 (misses sub-358m slips, keeps lower bound negative) DRAFT raw Tukey category_one_defects: upper=0.0 flags=202 (15.1%) DRAFT raw Tukey category_two_defects: upper=16.0 flags=43 (3.2%) DRAFT raw Tukey quakers: upper=0.0 flags=94 (7.0%) DRAFT MV(10 cols) flags n= 76 (cond # of cov: 67 ) NEW MV(7 cols) cond # of cov: 22 DRAFT abs(rstandard)>2 flags n= 57 (~ 4.3 % -- expected ~5% by construction)
command: create path: /app/outputs/audit_memo.md file_text: # Coffee Quality Outlier Audit , Methods Memo The draft applied one off-the-shelf recipe to every column. Each block below states why that recipe fails for the column's actual shape and what replaced it. All figures come from re-running `/app/analysis.R` on the 1,339-row CQI export. ## Data prep One withdrawn lot carries `0` for all ten grade components (and `total_cup_points = 0`). The draft's `total_cup_points > 0` filter happens to remove it, but keying on a derived total is fragile. I drop the row explicitly when **all ten** grade columns are zero , the documented sentinel , leaving `n_input = 1339` and `n_after_zero_score_drop = 1338`. Every downstream step runs on the 1,338 survivors in input order with a 0-indexed `row_id`. ## Altitude `altitude_mean_meters` is strongly right-skewed (median 1,311 m, max 190,164 m). A meter-scale Tukey fence with `k=3` returns `[-400, 3100]`: the lower bound is physically impossible, so it can never flag the many *low* decimal slips (`1`, `1.2`, `16.8` m), while its upper bound wrongly brands genuine highland lots (~3,280 m) as outliers. I instead build the fence on `log10(altitude)` over positive values (`k=3`) and back-transform, giving a symmetric, scale-appropriate `[357.4, 4923.8]` m that flags 51 rows. For each flagged row I test power-of-ten fixes on the first numeric token of the raw `altitude` string (`÷10`, `÷100`, then as-is) and keep the first candidate inside the fence; this recovers **7** displaced values (e.g. `190164`→1901.6, `11000 metros`→1100, `1100.00 mosl` whose stored mean was 110000→1100). Genuine low-altitude or unrecoverable rows keep `altitude_corrected_m = NA`. ## Defects `category_one_defects`, `category_two_defects`, and `quakers` are mass-at-zero (85%, 28%, 93% zeros) with a thin upper tail. Because `Q1=Q3=0` for two of them, the draft's raw IQR fence collapses to an upper bound of **0**, so it flags *every* nonzero lot: 202 rows (15.1%) for category one and 94 (7.0%) for quakers , far from "genuinely extreme." An IQR fence assumes a spread that a spike-at-zero distribution does not have. I use a per-column empirical **98th percentile** as the upper threshold and flag counts strictly above it, which isolates the true tail: thresholds 5 / 20 / 2 flag 20 / 23 / 25 rows (1.5% / 1.7% / 1.9%) , squarely the intended top 1–2%. A lot is a defect outlier if any column trips (63 rows). ## Multivariate grades The draft's Mahalanobis distance used all ten grade columns, including `uniformity`, `clean_cup`, and `sweetness`, which are near-constant (≈91% of lots score a perfect 10). Near-constant columns give the covariance matrix tiny eigenvalues, so their rare dips dominate the distance and destabilise the inverse (condition number ≈ 67). I restrict to the seven **primary SCA attributes** (`aroma`, `flavor`, `aftertaste`, `acidity`, `body`, `balance`, `cupper_points`), which are continuous and well-conditioned (condition number ≈ 22), and keep a classical mean/covariance Mahalanobis (as specified, not MCD). Against a `χ²(0.999, df=7) = 24.32` cut this flags 40 joint outliers , lots whose flavor *profile*, not just one score, is anomalous. ## Influence `abs(rstandard) > 2` measures residual size, not influence: by construction ~5% of any fit exceeds it (here 57 rows, 4.3%), and it ignores leverage entirely. The draft then copied that single count into leverage, Cook, and "both," which is simply wrong. I fit three OLS specs on the 1,338 grade-complete rows and pick the lowest finite **AIC**: grades_only (4459.4), grades_altitude (4461.3, `log10` elevation using corrected meters, mean-imputed), and grades_defects (**4418.7**, `log1p` defects, mean-imputed). grades_defects wins, confirming defect load carries real signal. On that fit (`n=1338`, `p=11`) I flag leverage `> 2p/n = 0.0164` (110 rows) and Cook's D `> 4/n = 0.0030` (66 rows) **separately**, and define influence as **both** (35 rows) , the standard …[truncated]
New file created successfully at /app/outputs/audit_memo.md
cd /app && sed -i 's/parametric Wald CID/parametric Wald CI/' outputs/audit_memo.md && echo "Word count: $(wc -w < outputs/audit_memo.md)" && echo "=== Headings present ===" && grep '^## ' outputs/audit_memo.md
Word count: 870 === Headings present === ## Data prep ## Altitude ## Defects ## Multivariate grades ## Influence ## Country ranking ## Sensitivity ## Composite flags
command: str_replace
path: /app/analysis.R
old_str: list(country = unbox(country), n = unbox(as.integer(n)),
raw_mean = unbox(round(raw_mean, 6)),
trimmed_mean = unbox(round(trimmed_mean, 6)),
new_str: list(country = unbox(country), n = unbox(as.integer(n)),
raw_mean = unbox(as.numeric(raw_mean)),
trimmed_mean = unbox(as.numeric(trimmed_mean)),The file /app/analysis.R has been edited successfully.
cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R && echo "--- validating schema ---" && Rscript -e '
suppressPackageStartupMessages(library(jsonlite))
j <- fromJSON("outputs/outlier_report.json", simplifyVector = FALSE)
req_top <- c("n_input","n_after_zero_score_drop","methods","altitude","defects","multivariate_grades","influence","country_ranking","sensitivity","flag_interactions")
cat("Top-level keys OK:", all(req_top %in% names(j)), "\n")
cat("methods keys OK:", all(c("altitude","defects","multivariate_grades","influence","country_ranking") %in% names(j$methods)), "\n")
cat("altitude keys OK:", all(c("lower_bound_m","upper_bound_m","n_flagged","n_unit_corrected") %in% names(j$altitude)), "\n")
cat("defects subkeys OK:", all(sapply(j$defects, function(d) all(c("upper_threshold","n_flagged") %in% names(d)))), "\n")
cat("mv keys OK:", all(c("columns_used","threshold","n_flagged") %in% names(j$multivariate_grades)), "\n")
cat("influence keys OK:", all(c("selected_model","aic_grades_only","aic_grades_altitude","aic_grades_defects","leverage_threshold","cook_threshold","n_high_leverage","n_high_cook","n_both","top10_cook_row_ids") %in% names(j$influence)), "\n")
cat("country keys OK:", all(c("n_countries_evaluated","n_rank_changes_ge_2","max_abs_rank_change","per_country") %in% names(j$country_ranking)), "\n")
pc1 <- j$country_ranking$per_country[[1]]
cat("per_country item keys OK:", all(c("country","n","raw_mean","trimmed_mean","raw_rank","trimmed_rank") %in% names(pc1)), "\n")
cat("sensitivity keys OK:", all(c("coef_flavor_full","coef_flavor_dropped","delta_pct","delta_pct_ci_lower","delta_pct_ci_upper","n_bootstrap_replicates") %in% names(j$sensitivity)), "\n")
cat("flag_interactions keys OK:", all(c("n_with_zero_flags","n_with_one_flag","n_with_two_flags","n_with_three_flags","n_with_four_flags","n_alt_only","n_mv_only","n_influence_only","n_defect_only","n_alt_and_mv","n_alt_and_influence","n_mv_and_influence","n_alt_mv_influence") %in% names(j$flag_interactions)), "\n")
cat("selected_model valid:", j$influence$selected_model %in% c("grades_only","grades_altitude","grades_defects"), "\n")
cat("top10 length:", length(j$influence$top10_cook_row_ids), " columns_used length:", length(j$multivariate_grades$columns_used), "\n")
cat("per_country count:", length(j$country_ranking$per_country), "\n")
cat("Ethiopia raw_mean full precision:", format(j$country_ranking$per_country[[1]]$raw_mean, digits=15), "\n")
'Audit complete: n_input=1339, n_after=1338, selected=grades_defects, composite=105 --- validating schema --- Top-level keys OK: TRUE methods keys OK: TRUE altitude keys OK: TRUE defects subkeys OK: TRUE mv keys OK: TRUE influence keys OK: TRUE country keys OK: TRUE per_country item keys OK: TRUE sensitivity keys OK: TRUE flag_interactions keys OK: TRUE selected_model valid: TRUE top10 length: 10 columns_used length: 7 per_country count: 21 Ethiopia raw_mean full precision: 85.4840909091
cd /app && python3 -c "
import json
j = json.load(open('outputs/outlier_report.json'))
def check(path, v):
import numbers
return isinstance(v, (numbers.Number, str, bool))
# scalars must be bare
scalars = [j['n_input'], j['n_after_zero_score_drop'], j['altitude']['lower_bound_m'],
j['influence']['selected_model'], j['influence']['aic_grades_defects'],
j['sensitivity']['delta_pct'], j['flag_interactions']['n_alt_only']]
print('all scalars bare (not list):', all(not isinstance(x, list) for x in scalars))
print('columns_used is list:', isinstance(j['multivariate_grades']['columns_used'], list))
print('top10 is list of len 10:', isinstance(j['influence']['top10_cook_row_ids'], list), len(j['influence']['top10_cook_row_ids']))
print('per_country is list:', isinstance(j['country_ranking']['per_country'], list))
print('n_input =', j['n_input'], '| selected =', j['influence']['selected_model'])
print('aic fields:', j['influence']['aic_grades_only'], j['influence']['aic_grades_altitude'], j['influence']['aic_grades_defects'])
print('JSON parses cleanly, no NaN/Inf tokens:', 'NaN' not in open('outputs/outlier_report.json').read())
"all scalars bare (not list): True columns_used is list: True top10 is list of len 10: True 10 per_country is list: True n_input = 1339 | selected = grades_defects aic fields: 4459.3615159695 4461.292329589 4418.6800742954 JSON parses cleanly, no NaN/Inf tokens: True
cd /app && python3 -c "
import pandas as pd
df = pd.read_csv('outputs/outlier_flags.csv')
print('shape:', df.shape)
print('dtypes:'); print(df.dtypes)
print()
print('row_id range:', df.row_id.min(), '->', df.row_id.max(), '| monotonic:', df.row_id.is_monotonic_increasing)
print('altitude_corrected_m non-null count:', df.altitude_corrected_m.notna().sum())
print('flag sums: alt=%d defect=%d mv=%d lev=%d cook=%d composite=%d' % (
df.altitude_outlier_flag.sum(), df.defect_outlier_flag.sum(),
df.multivariate_grade_outlier_flag.sum(), df.high_leverage_flag.sum(),
df.high_cook_flag.sum(), df.composite_drop_flag.sum()))
print('bool dtype flags:', all(df[c].dtype==bool for c in ['altitude_outlier_flag','defect_outlier_flag','multivariate_grade_outlier_flag','high_leverage_flag','high_cook_flag','composite_drop_flag']))
# composite consistency in pandas
comp = df.altitude_outlier_flag | df.multivariate_grade_outlier_flag | (df.high_leverage_flag & df.high_cook_flag)
print('composite matches:', (comp == df.composite_drop_flag).all())
"shape: (1338, 10) dtypes: row_id int64 country_of_origin object total_cup_points float64 altitude_outlier_flag bool altitude_corrected_m float64 defect_outlier_flag bool multivariate_grade_outlier_flag bool high_leverage_flag bool high_cook_flag bool composite_drop_flag bool dtype: object row_id range: 0 -> 1337 | monotonic: True altitude_corrected_m non-null count: 7 flag sums: alt=51 defect=63 mv=40 lev=110 cook=66 composite=105 bool dtype flags: True composite matches: True
cd /tmp && Rscript -e '
# Fresh session: only the rds, nothing else from the analysis
fit <- readRDS("/app/outputs/coffee_ols_selected.rds")
n <- nobs(fit); p <- length(coef(fit))
lev <- hatvalues(fit); ck <- cooks.distance(fit)
cat("Standalone reload OK. n =", n, "p =", p, "\n")
cat("n_high_lev =", sum(lev > 2*p/n), " n_high_cook =", sum(ck > 4/n), " n_both =", sum(lev>2*p/n & ck>4/n), "\n")
cat("coef names:", paste(names(coef(fit)), collapse=", "), "\n")
cat("flavor coef in selected fit:", round(coef(fit)["flavor"],4), "\n")
'Standalone reload OK. n = 1338 p = 11 n_high_lev = 110 n_high_cook = 66 n_both = 35 coef names: (Intercept), aroma, flavor, aftertaste, acidity, body, balance, cupper_points, l1p_category_one_defects, l1p_category_two_defects, l1p_quakers flavor coef in selected fit: 2.0801
cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R && echo "=== Final outputs ===" && ls -la outputs/
Audit complete: n_input=1339, n_after=1338, selected=grades_defects, composite=105 === Final outputs === total 226 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:26 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:25 .. -rw-r--r-- 1 h2tagent h2tagent 137312 Jun 20 18:26 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78269 Jun 20 18:26 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6655 Jun 20 18:26 outlier_report.json
command: str_replace
path: /app/analysis.R
old_str: write_json(results, file.path(OUTPUT_DIR, "outlier_report.json"),
auto_unbox = FALSE, pretty = TRUE, na = "null", digits = 10)
new_str: write_json(results, file.path(OUTPUT_DIR, "outlier_report.json"),
auto_unbox = FALSE, pretty = TRUE, na = "null", digits = 10)
# ---------------------------------------------------------------------------
# audit_memo.md -- generated so a clean re-run reproduces this artifact too.
# Draft-comparison figures are recomputed here to keep the prose faithful.
# ---------------------------------------------------------------------------
draft_q <- quantile(alt_pos, c(0.25, 0.75))
draft_lo <- draft_q[[1]] - 3 * (draft_q[[2]] - draft_q[[1]])
draft_hi <- draft_q[[2]] + 3 * (draft_q[[2]] - draft_q[[1]])
draft_defect_tukey <- function(x) {
q <- quantile(x, c(0.25, 0.75), na.rm = TRUE)
up <- q[[2]] + 3 * (q[[2]] - q[[1]]); sum(!is.na(x) & x > up)
}
d1_draft <- draft_defect_tukey(dat$category_one_defects)
dq_draft <- draft_defect_tukey(dat$quakers)
G10 <- as.matrix(dat[, GRADES10]); cc10 <- complete.cases(G10)
S10 <- cov(G10[cc10, ]); md10 <- mahalanobis(G10[cc10, ], colMeans(G10[cc10, ]), S10, tol = 1e-30)
mv10_draft <- sum(md10 > qchisq(0.999, length(GRADES10)))
kappa10 <- round(kappa(S10)); kappa7 <- round(kappa(mv_S))
rstd_draft <- sum(abs(rstandard(m_go)) > 2)
pc1 <- country_tbl[1, ]
movers <- country_tbl$country[country_tbl$abs_rank_change >= 2]
nf <- function(x, d = 1) formatdec(x, d)
formatdec <- function(x, d = 1) formatC(x, format = "f", digits = d, big.mark = "")
f0 <- function(x) formatC(round(x), format = "d", big.mark = "")
memo <- paste0(
"# Coffee Quality Outlier Audit \u2014 Methods Memo\n\n",
"The draft applied one off-the-shelf recipe to every column. Each block below ",
"states why that recipe fails for the column's actual shape and what replaced ",
"it. Figures are regenerated by `/app/analysis.R`.\n\n",
"## Data prep\n",
"One withdrawn lot records `0` for all ten grade components (and ",
"`total_cup_points = 0`). The draft's `total_cup_points > 0` filter happens to ",
"remove it, but keying on a derived total is fragile. I drop the row explicitly ",
"when **all ten** grade columns are zero \u2014 the documented sentinel \u2014 leaving ",
"`n_input = ", f0(n_input), "` and `n_after_zero_score_drop = ", f0(n_after),
"`. Everything downstream runs on the ", f0(n_after), " survivors in input ",
"order with a 0-indexed `row_id`.\n\n",
"## Altitude\n",
"`altitude_mean_meters` is strongly right-skewed (median ~1,311 m, max ",
"190,164 m). A meter-scale Tukey fence (k=3) returns `[", f0(draft_lo), ", ",
f0(draft_hi), "]`: the lower bound is physically impossible, so it can never ",
"flag the many *low* decimal slips (`1`, `1.2`, `16.8` m), while its upper ",
"bound brands genuine highland lots as outliers. I build the fence on ",
"`log10(altitude)` over positive values (k=3) and back-transform, giving a ",
"scale-appropriate `[", nf(alt_lo_m), ", ", nf(alt_hi_m), "]` m that flags ",
f0(n_flag_alt), " rows. For each flagged row I test power-of-ten fixes on the ",
"first numeric token of the raw string (\u00f710, \u00f7100, then as-is) and keep the ",
"first candidate inside the fence; this recovers **", f0(n_unit_corrected),
"** displaced values (e.g. `190164`\u21921901.6, `11000 metros`\u21921100, ",
"`1100.00 mosl` whose stored mean was 110000\u21921100). Unrecoverable or ",
"genuinely low rows keep `altitude_corrected_m = NA`.\n\n",
"## Defects\n",
"`category_one_defects`, `category_two_defects` and `quakers` are mass-at-zero ",
"(85%, 28%, 93% zeros) with a thin upper tail. Because `Q1=Q3=0` for two of ",
"them, the draft's raw IQR fence collapses to an upper bound of **0** and flags ",
"*every* nonzero lot: ", f0(d1_draft), " rows for category one and ", f0(dq_draft),
" for quakers \u2014 nowhere near 'genuinely extreme'. An IQR fence assumes a spread ",
"a spike-at-zero distribution lacks. I use a per-c
…[truncated]The file /app/analysis.R has been edited successfully.
command: str_replace path: /app/analysis.R old_str: nf <- function(x, d = 1) formatdec(x, d) formatdec <- function(x, d = 1) formatC(x, format = "f", digits = d, big.mark = "") f0 <- function(x) formatC(round(x), format = "d", big.mark = "") new_str: nf <- function(x, d = 1) formatC(x, format = "f", digits = d, big.mark = "") f0 <- function(x) formatC(round(x), format = "d", big.mark = "")
The file /app/analysis.R has been edited successfully.
cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R && echo "=== outputs ===" && ls -la outputs/ && echo "=== memo word count ===" && wc -w outputs/audit_memo.md && echo "=== headings ===" && grep '^## ' outputs/audit_memo.md
Audit complete: n_input=1339, n_after=1338, selected=grades_defects, composite=105 === outputs === total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:26 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:26 .. -rw-r--r-- 1 h2tagent h2tagent 5115 Jun 20 18:26 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137312 Jun 20 18:26 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78269 Jun 20 18:26 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6655 Jun 20 18:26 outlier_report.json === memo word count === 781 outputs/audit_memo.md === headings === ## Data prep ## Altitude ## Defects ## Multivariate grades ## Influence ## Country ranking ## Sensitivity ## Composite flags
cd /app && cat outputs/audit_memo.md
# Coffee Quality Outlier Audit , Methods Memo The draft applied one off-the-shelf recipe to every column. Each block below states why that recipe fails for the column's actual shape and what replaced it. Figures are regenerated by `/app/analysis.R`. ## Data prep One withdrawn lot records `0` for all ten grade components (and `total_cup_points = 0`). The draft's `total_cup_points > 0` filter happens to remove it, but keying on a derived total is fragile. I drop the row explicitly when **all ten** grade columns are zero , the documented sentinel , leaving `n_input = 1339` and `n_after_zero_score_drop = 1338`. Everything downstream runs on the 1338 survivors in input order with a 0-indexed `row_id`. ## Altitude `altitude_mean_meters` is strongly right-skewed (median ~1,311 m, max 190,164 m). A meter-scale Tukey fence (k=3) returns `[-400, 3100]`: the lower bound is physically impossible, so it can never flag the many *low* decimal slips (`1`, `1.2`, `16.8` m), while its upper bound brands genuine highland lots as outliers. I build the fence on `log10(altitude)` over positive values (k=3) and back-transform, giving a scale-appropriate `[357.4, 4923.8]` m that flags 51 rows. For each flagged row I test power-of-ten fixes on the first numeric token of the raw string (÷10, ÷100, then as-is) and keep the first candidate inside the fence; this recovers **7** displaced values (e.g. `190164`→1901.6, `11000 metros`→1100, `1100.00 mosl` whose stored mean was 110000→1100). Unrecoverable or genuinely low rows keep `altitude_corrected_m = NA`. ## Defects `category_one_defects`, `category_two_defects` and `quakers` are mass-at-zero (85%, 28%, 93% zeros) with a thin upper tail. Because `Q1=Q3=0` for two of them, the draft's raw IQR fence collapses to an upper bound of **0** and flags *every* nonzero lot: 202 rows for category one and 94 for quakers , nowhere near 'genuinely extreme'. An IQR fence assumes a spread a spike-at-zero distribution lacks. I use a per-column empirical **98th percentile** as the upper cut and flag counts strictly above it: thresholds 5 / 20 / 2 flag 20 / 23 / 25 rows (~1.5–1.9%) , the intended top 1–2%. A lot is a defect outlier if any column trips. ## Multivariate grades The draft's Mahalanobis used all ten grades, including `uniformity`, `clean_cup` and `sweetness`, which are near-constant (~91% of lots score a perfect 10). Near-constant columns give the covariance tiny eigenvalues, so their rare dips dominate the distance and destabilise the inverse (condition number ~67), flagging 76 rows. I restrict to the seven **primary SCA attributes**, which are continuous and far better conditioned (~22), and keep a classical mean/covariance Mahalanobis (not MCD). Against `χ²(0.999, df=7) = 24.32` this flags 40 lots whose flavor *profile*, not one score, is anomalous. ## Influence `abs(rstandard) > 2` measures residual size, not influence: ~5% of any fit exceeds it by construction (here 57 rows) and it ignores leverage entirely. The draft then copied that one count into leverage, Cook and 'both', which is wrong. I fit three OLS specs on the 1338 grade-complete rows and pick the lowest finite **AIC**: grades_only (4459.4), grades_altitude (4461.3), grades_defects (**4418.7**). grades_defects wins, confirming defect load carries real signal. On it (n=1338, p=11) I flag leverage `> 2p/n = 0.0164` (110 rows) and Cook's D `> 4/n = 0.0030` (66 rows) **separately**, and define influence as **both** (35 rows). ## Country ranking Raw country means let a few extreme lots swing small samples, and the draft reported no robustness check (trimmed mean = raw mean). For the 21 countries with ≥10 lots I compare the raw mean with a **20% trimmed mean** of `total_cup_points` and rank both. 2 countries shift ≥2 positions (Taiwan, Guatemala), max absolute move 2 , small, but it exposes which rankings depend on tail lots rather than typical quality. ## Sensitivity The draft never quantified how much outliers move a coefficient. I define the …[truncated]
cd /app && mkdir -p /tmp/run_a && cp outputs/* /tmp/run_a/ && rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1 && echo "Comparing two independent clean runs:" && for f in outlier_report.json outlier_flags.csv audit_memo.md; do if diff -q /tmp/run_a/$f outputs/$f >/dev/null; then echo " $f: IDENTICAL"; else echo " $f: DIFFERS"; fi; done && echo "RDS coef check:" && Rscript -e 'a<-readRDS("/tmp/run_a/coffee_ols_selected.rds"); b<-readRDS("outputs/coffee_ols_selected.rds"); cat(" coefs identical:", isTRUE(all.equal(coef(a),coef(b))), "\n")'Comparing two independent clean runs: outlier_report.json: IDENTICAL outlier_flags.csv: IDENTICAL audit_memo.md: IDENTICAL RDS coef check: coefs identical: TRUE
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
cat("raw rows:", nrow(raw), " csv rows:", nrow(csv), " dropped:", nrow(raw)-nrow(csv), "\n")
cat("NA-country rows retained in CSV:", sum(is.na(csv$country_of_origin)), "\n")
cat("Any total_cup_points==0 remaining in CSV:", sum(csv$total_cup_points==0), "\n")
cat("min total_cup_points in CSV:", min(csv$total_cup_points), "\n")
# Confirm dropped row was the all-zero sentinel
g10 <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
dropped <- raw[rowSums(raw[,g10]==0, na.rm=TRUE)==10, ]
cat("Sentinel dropped had total_cup_points:", dropped$total_cup_points, "owner:", dropped$owner, "\n")
' 2>&1 | grep -v "^Warning\|parsing\|Column\|expected\|actual\|file\|Specify\|dat <-\|Rows:\|•\|──\|ℹ\|See problems\|Delim\|chr\|dbl\|date"raw rows: 1339 csv rows: 1338 dropped: 1 NA-country rows retained in CSV: 1 Any total_cup_points==0 remaining in CSV: 0 min total_cup_points in CSV: 59.83 Sentinel dropped had total_cup_points: 0 owner: bismarck castro
cd /app && sed -n '1,130p' analysis.R
#!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-matched rewrite
# =====================================================
# Each column is audited with a method that matches its distributional shape,
# replacing the draft's off-the-shelf summaries. See audit_memo.md for the
# rationale behind every rejection. This script alone reproduces all artifacts
# from a clean /app/outputs/ directory.
suppressPackageStartupMessages({
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(jsonlite)
library(purrr)
})
DATA_PATH <- "/app/data/coffee_ratings.csv"
OUTPUT_DIR <- "/app/outputs"
dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE)
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
# ---------------------------------------------------------------------------
# 1. Sentinel drop -- one withdrawn lot has every grade recorded as 0.
# ---------------------------------------------------------------------------
GRADES10 <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"uniformity", "clean_cup", "sweetness", "cupper_points")
all_zero <- rowSums(raw[, GRADES10] == 0, na.rm = TRUE) == length(GRADES10)
dat <- raw[!all_zero, , drop = FALSE]
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L))
N <- nrow(dat)
# ---------------------------------------------------------------------------
# 2. Altitude -- log10-scale Tukey fence (k = 3); meter-scale fences mislead
# on this right-skewed column. Recover decimal-displacement unit slips.
# ---------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
Lq <- quantile(log10(alt_pos), c(0.25, 0.75))
L_iqr <- Lq[[2]] - Lq[[1]]
alt_lo_m <- 10^(Lq[[1]] - 3 * L_iqr)
alt_hi_m <- 10^(Lq[[2]] + 3 * L_iqr)
flag_alt <- !is.na(alt) & alt > 0 & (alt < alt_lo_m | alt > alt_hi_m)
# First numeric token of the raw altitude string (digits + optional decimals;
# thousands-separators/units are intentionally not treated as part of it).
first_token <- function(s) as.numeric(str_extract(s, "[0-9]+\\.?[0-9]*"))
alt_corrected <- rep(NA_real_, N)
for (i in which(flag_alt)) {
tok <- first_token(dat$altitude[i])
if (is.na(tok)) next
for (cand in c(tok / 10, tok / 100, tok)) { # plausible power-of-ten fixes
if (!is.na(cand) && cand >= alt_lo_m && cand <= alt_hi_m) {
alt_corrected[i] <- cand
break
}
}
}
n_flag_alt <- sum(flag_alt)
n_unit_corrected <- sum(!is.na(alt_corrected))
# ---------------------------------------------------------------------------
# 3. Defects -- mass-at-zero counts with a sparse upper tail. Raw Tukey fences
# collapse to 0 (Q1 = Q3 = 0) and flag every nonzero row, so use a per-column
# empirical 98th-percentile upper threshold; flag counts strictly above it.
# ---------------------------------------------------------------------------
DEFECT_COLS <- c("category_one_defects", "category_two_defects", "quakers")
defect_flag_one <- function(x) {
thr <- unname(quantile(x, 0.98, na.rm = TRUE, type = 7))
list(upper_threshold = thr,
flag_vec = !is.na(x) & x > thr,
n = sum(!is.na(x) & x > thr))
}
defect_res <- lapply(DEFECT_COLS, function(cn) defect_flag_one(dat[[cn]]))
names(defect_res) <- DEFECT_COLS
flag_defect <- Reduce(`|`, lapply(defect_res, `[[`, "flag_vec"))
# ---------------------------------------------------------------------------
# 4. Multivariate grades -- classical Mahalanobis on the SEVEN primary SCA
# attributes only (the cleanliness columns are near-constant and inflate /
# destabilise the covariance). Chi-square upper tail at 99.9%, df = 7.
# ---------------------------------------------------------------------------
MV_COLS <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"cupper_points")
G <- as.matrix(dat[, MV_COLS])
G_cc <- complete.ca
…[truncated]cd /app && sed -n '130,230p' analysis.R
m_gd <- lm(as.formula(paste("total_cup_points ~", GFORM, "+",
paste(defx_names, collapse = " + "))), data = fr)
aic_vals <- c(grades_only = AIC(m_go),
grades_altitude = AIC(m_ga),
grades_defects = AIC(m_gd))
finite_aic <- aic_vals[is.finite(aic_vals)]
selected_model <- if (length(finite_aic) == 0) "grades_only" else
names(finite_aic)[which.min(finite_aic)]
sel_fit <- switch(selected_model,
grades_only = m_go, grades_altitude = m_ga, grades_defects = m_gd)
n_fit <- nobs(sel_fit)
p_fit <- length(coef(sel_fit))
lev_thr <- 2 * p_fit / n_fit
cook_thr <- 4 / n_fit
lev_vec <- hatvalues(sel_fit)
cook_vec <- cooks.distance(sel_fit)
hi_lev_local <- lev_vec > lev_thr
hi_cook_local <- cook_vec > cook_thr
both_local <- hi_lev_local & hi_cook_local
# expand fit-frame diagnostics back to full post-drop row space
hi_lev <- rep(FALSE, N); hi_lev[reg_idx] <- hi_lev_local
hi_cook <- rep(FALSE, N); hi_cook[reg_idx] <- hi_cook_local
flag_infl <- rep(FALSE, N); flag_infl[reg_idx] <- both_local
# ten largest Cook's D -> 0-indexed post-drop row ids (selected fit)
top_local <- order(cook_vec, decreasing = TRUE)[1:10]
top10_row_ids <- as.integer(dat$row_id[reg_idx[top_local]])
# ---------------------------------------------------------------------------
# 6. Country ranking -- raw mean vs 20% trimmed (robust) mean, >= 10 lots.
# ---------------------------------------------------------------------------
country_tbl <- dat %>%
filter(!is.na(country_of_origin)) %>%
group_by(country = country_of_origin) %>%
summarise(n = n(),
raw_mean = mean(total_cup_points),
trimmed_mean = mean(total_cup_points, trim = 0.2),
.groups = "drop") %>%
filter(n >= 10) %>%
mutate(raw_rank = rank(-raw_mean, ties.method = "min"),
trimmed_rank = rank(-trimmed_mean, ties.method = "min"),
abs_rank_change = abs(raw_rank - trimmed_rank)) %>%
arrange(raw_rank, country)
n_countries <- nrow(country_tbl)
n_changes_ge_2 <- sum(country_tbl$abs_rank_change >= 2)
max_change <- if (n_countries > 0) max(country_tbl$abs_rank_change) else 0L
# ---------------------------------------------------------------------------
# 7. Sensitivity -- composite = altitude OR multivariate OR (leverage AND cook).
# Defect outliers are deliberately excluded from the composite.
# ---------------------------------------------------------------------------
composite <- flag_alt | flag_mv | flag_infl
# refit grades_only on regression frame before/after dropping composite set
sens_full_df <- dat[reg_idx, , drop = FALSE]
sens_keep <- !composite[reg_idx]
m_sens_full <- lm(f_go, data = sens_full_df)
m_sens_drop <- lm(f_go, data = sens_full_df[sens_keep, , drop = FALSE])
coef_full <- unname(coef(m_sens_full)["flavor"])
coef_dropped <- unname(coef(m_sens_drop)["flavor"])
delta_pct <- 100 * (coef_dropped - coef_full) / coef_full
# percentile bootstrap CI on the percent change (carry ORIGINAL composite flags)
set.seed(20260512)
B <- 500
boot_n <- nrow(sens_full_df)
comp_frame <- composite[reg_idx]
boot_delta <- numeric(0)
for (b in seq_len(B)) {
idx <- sample.int(boot_n, boot_n, replace = TRUE)
d_b <- sens_full_df[idx, , drop = FALSE]
comp_b <- comp_frame[idx] # original flags, not recomputed
keep_b <- d_b[!comp_b, , drop = FALSE]
if (nrow(keep_b) < 10) next
ff <- tryCatch(lm(f_go, data = d_b), error = function(e) NULL)
fd <- tryCatch(lm(f_go, data = keep_b), error = function(e) NULL)
if (is.null(ff) || is.null(fd)) next
cf <- unname(coef(ff)["flavor"]); cd <- unname(coef(fd)["flavor"])
if (is.na(cf) || is.na(cd) || cf == 0) next
boot_delta <- c(boot_delta, 100 * (cd - cf) / cf)
}
n_boot <- length(boot_delta)
if (n_boot > 0) {
ci <- unname(quantile(boot_delta, c(0.025, 0.975)))
ci_lower <- ci[1]; ci_upper <- ci[2]
} else {
ci_lower <- NA_real_; ci_upper <- NA_real_
}
#
…[truncated]cd /app && rm -rf outputs && Rscript analysis.R && echo "EXIT CODE: $?" && echo "=== All invariants ===" && Rscript -e '
suppressPackageStartupMessages({library(jsonlite); library(readr)})
j <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
fit <- readRDS("outputs/coffee_ols_selected.rds")
ok <- function(label, cond) cat(sprintf("[%s] %s\n", ifelse(cond,"PASS","FAIL"), label))
ok("n_input=1339", j$n_input==1339)
ok("n_after=1338", j$n_after_zero_score_drop==1338)
ok("csv rows=1338", nrow(csv)==1338)
ok("alt bounds ordered & positive", j$altitude$lower_bound_m>0 && j$altitude$lower_bound_m < j$altitude$upper_bound_m)
ok("alt n_flagged matches csv", j$altitude$n_flagged==sum(csv$altitude_outlier_flag))
ok("alt corrected matches csv non-NA", j$altitude$n_unit_corrected==sum(!is.na(csv$altitude_corrected_m)))
ok("corrected rows all altitude-flagged", all(csv$altitude_outlier_flag[!is.na(csv$altitude_corrected_m)]))
ok("corrected values within fence", all(csv$altitude_corrected_m[!is.na(csv$altitude_corrected_m)] >= j$altitude$lower_bound_m & csv$altitude_corrected_m[!is.na(csv$altitude_corrected_m)] <= j$altitude$upper_bound_m))
for(d in c("category_one_defects","category_two_defects","quakers")){
r <- 100*j$defects[[d]]$n_flagged/1338
ok(sprintf("defect %s in 1-2%% (%.2f%%)",d,r), r>=1 && r<=2)
}
ok("defect any matches csv", sum(csv$defect_outlier_flag)==63)
ok("mv 7 cols", length(j$multivariate_grades$columns_used)==7 && !("clean_cup" %in% unlist(j$multivariate_grades$columns_used)))
ok("mv threshold=qchisq(.999,7)", abs(j$multivariate_grades$threshold - qchisq(.999,7))<1e-9)
ok("mv n matches csv", j$multivariate_grades$n_flagged==sum(csv$multivariate_grade_outlier_flag))
ok("selected=grades_defects (lowest AIC)", j$influence$selected_model=="grades_defects")
ok("AIC selected is min finite", which.min(c(j$influence$aic_grades_only,j$influence$aic_grades_altitude,j$influence$aic_grades_defects))==3)
n<-nobs(fit);p<-length(coef(fit))
ok("lev_thr=2p/n", abs(j$influence$leverage_threshold-2*p/n)<1e-12)
ok("cook_thr=4/n", abs(j$influence$cook_threshold-4/n)<1e-12)
ok("lev flag recompute matches csv", all((hatvalues(fit)>2*p/n)==csv$high_leverage_flag))
ok("cook flag recompute matches csv", all((cooks.distance(fit)>4/n)==csv$high_cook_flag))
ok("n_both = lev&cook", j$influence$n_both==sum(csv$high_leverage_flag & csv$high_cook_flag))
ok("top10 cook correct", all(unlist(j$influence$top10_cook_row_ids)==csv$row_id[order(-cooks.distance(fit))[1:10]]))
ok("composite=alt|mv|(lev&cook)", all(csv$composite_drop_flag==(csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag))))
ok("defects NOT in composite (some defect-only exist)", any(csv$defect_outlier_flag & !csv$composite_drop_flag))
ok("country n_eval=21", j$country_ranking$n_countries_evaluated==21)
ok("delta_pct=100*(drop-full)/full", abs(j$sensitivity$delta_pct - 100*(j$sensitivity$coef_flavor_dropped-j$sensitivity$coef_flavor_full)/j$sensitivity$coef_flavor_full)<1e-6)
ok("CI brackets & ordered", j$sensitivity$delta_pct_ci_lower < j$sensitivity$delta_pct_ci_upper)
ok("bootstrap reps=500", j$sensitivity$n_bootstrap_replicates==500)
fi<-j$flag_interactions
ok("interaction buckets sum=1338", fi$n_with_zero_flags+fi$n_with_one_flag+fi$n_with_two_flags+fi$n_with_three_flags+fi$n_with_four_flags==1338)
' 2>&1 | grep -v "Rows:\|Columns\|Delim\|chr\|dbl\|lgl\|──\|ℹ\|•"Audit complete: n_input=1339, n_after=1338, selected=grades_defects, composite=105 EXIT CODE: 0 === All invariants === [PASS] n_input=1339 [PASS] n_after=1338 [PASS] csv rows=1338 [PASS] alt bounds ordered & positive [PASS] alt n_flagged matches csv [PASS] alt corrected matches csv non-NA [PASS] corrected rows all altitude-flagged [PASS] corrected values within fence [PASS] defect category_one_defects in 1-2% (1.49%) [PASS] defect category_two_defects in 1-2% (1.72%) [PASS] defect quakers in 1-2% (1.87%) [PASS] defect any matches csv [PASS] mv 7 cols [PASS] mv threshold=qchisq(.999,7) [PASS] mv n matches csv [PASS] selected=grades_defects (lowest AIC) [PASS] AIC selected is min finite [PASS] lev_thr=2p/n [PASS] cook_thr=4/n [PASS] lev flag recompute matches csv [PASS] cook flag recompute matches csv [PASS] n_both = lev&cook [PASS] top10 cook correct [PASS] composite=alt|mv|(lev&cook) [PASS] defects NOT in composite (some defect-only exist) [PASS] country n_eval=21 [PASS] delta_pct=100*(drop-full)/full [PASS] CI brackets & ordered [PASS] bootstrap reps=500 [PASS] interaction buckets sum=1338
rm -rf /tmp/run_a /tmp/run1.json /tmp/run1.csv /tmp/*.rds 2>/dev/null; cd /app && echo "=== Final /app/outputs/ contents ===" && ls -la outputs/ && echo && echo "=== Final sanity: valid JSON + all 4 artifacts present ===" && python3 -c "import json; json.load(open('outputs/outlier_report.json')); print('JSON valid')" && for f in outlier_report.json coffee_ols_selected.rds outlier_flags.csv audit_memo.md; do test -s outputs/$f && echo "present & non-empty: $f"; done=== Final /app/outputs/ contents === total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:27 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:27 .. -rw-r--r-- 1 h2tagent h2tagent 5115 Jun 20 18:27 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137312 Jun 20 18:27 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78269 Jun 20 18:27 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6655 Jun 20 18:27 outlier_report.json === Final sanity: valid JSON + all 4 artifacts present === JSON valid present & non-empty: outlier_report.json present & non-empty: coffee_ols_selected.rds present & non-empty: outlier_flags.csv present & non-empty: audit_memo.md
No file changes, the agent left the workspace untouched.
# SyncVals normalized verifier summary # task: coffee-ratings-outliers # attempt: 8 status: completed raw_score: 0.0 # reward rule: binary: score 0.0 != 1.0 # => reward = 0.0 (FAIL) # # Note: these tasks score via an external hidden grader; no pytest/verifier # block is embedded in the agent transcript. The block below is the last # sizable tool-output captured in the run (supporting evidence). ---------------------------------------------------------------------- === Final /app/outputs/ contents === total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:27 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:27 .. -rw-r--r-- 1 h2tagent h2tagent 5115 Jun 20 18:27 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137312 Jun 20 18:27 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78269 Jun 20 18:27 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6655 Jun 20 18:27 outlier_report.json === Final sanity: valid JSON + all 4 artifacts present === JSON valid present & non-empty: outlier_report.json present & non-empty: coffee_ols_selected.rds present & non-empty: outlier_flags.csv present & non-empty: audit_memo.md
Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_026eff7e92414f77. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_026eff7e92414f77 · verifier authoritative; classifier explanatory.