tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.cd /app && ls -la && echo "---OUTPUTS---" && 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--- 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 -3 data/coffee_ratings.csv && echo "---COLS---" && head -1 data/coffee_ratings.csv | tr ',' '\n' | nl
total_cup_points,species,owner,country_of_origin,farm_name,lot_number,mill,ico_number,company,altitude,region,producer,number_of_bags,bag_weight,in_country_partner,harvest_year,grading_date,owner_1,variety,processing_method,aroma,flavor,aftertaste,acidity,body,balance,uniformity,clean_cup,sweetness,cupper_points,moisture,category_one_defects,quakers,color,category_two_defects,expiration,certification_body,certification_address,certification_contact,unit_of_measurement,altitude_low_meters,altitude_high_meters,altitude_mean_meters
90.58,Arabica,metad plc,Ethiopia,metad plc,NA,metad plc,2014/2015,metad agricultural developmet plc,1950-2200,guji-hambela,METAD PLC,300,60 kg,METAD Agricultural Development plc,2014,"April 4th, 2015",metad plc,NA,Washed / Wet,8.67,8.83,8.67,8.75,8.5,8.42,10,10,10,8.75,0.12,0,0,Green,0,"April 3rd, 2016",METAD Agricultural Development plc,309fcf77415a3661ae83e027f7e5f05dad786e44,19fef5a731de2db57d16da10287413f5f99bc2dd,m,1950,2200,2075
89.92,Arabica,metad plc,Ethiopia,metad plc,NA,metad plc,2014/2015,metad agricultural developmet plc,1950-2200,guji-hambela,METAD PLC,300,60 kg,METAD Agricultural Development plc,2014,"April 4th, 2015",metad plc,Other,Washed / Wet,8.75,8.67,8.5,8.58,8.42,8.42,10,10,10,8.58,0.12,0,0,Green,1,"April 3rd, 2016",METAD Agricultural Development plc,309fcf77415a3661ae83e027f7e5f05dad786e44,19fef5a731de2db57d16da10287413f5f99bc2dd,m,1950,2200,2075
---COLS---
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 && which Rscript && Rscript --version && echo "---PKGS---" && Rscript -e 'for (p in c("dplyr","tidyr","readr","stringr","jsonlite","purrr","MASS")) cat(p, as.character(packageVersion(p)), "\n")'/usr/local/bin/Rscript Rscript (R) version 4.4.2 (2024-10-31) ---PKGS--- dplyr 1.1.4 tidyr 1.3.1 readr 2.1.5 stringr 1.5.1 jsonlite 1.9.0 purrr 1.0.4 MASS 7.3.64
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","total_cup_points")
# zero-score lots
z <- raw %>% filter(aroma==0 & flavor==0 & aftertaste==0 & acidity==0 & body==0 & balance==0 & cupper_points==0)
cat("rows with all 7 primary grades zero:", nrow(z), "\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[raw$total_cup_points==0, c("total_cup_points","aroma","flavor","acidity","country_of_origin")])
'n_input = 1339
rows with all 7 primary grades zero: 1
rows with total_cup_points==0: 1
rows with total_cup_points<=0: 1
# A tibble: 1 × 5
total_cup_points aroma flavor acidity country_of_origin
<dbl> <dbl> <dbl> <dbl> <chr>
1 0 0 0 0 Hondurascd /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)
cat("n_after =", nrow(dat), "\n\n")
cat("=== altitude_mean_meters summary ===\n")
print(summary(dat$altitude_mean_meters))
cat("NA count:", sum(is.na(dat$altitude_mean_meters)), "\n")
cat("zero/neg count:", sum(dat$altitude_mean_meters<=0, na.rm=TRUE), "\n\n")
a <- dat$altitude_mean_meters
cat("top values:\n"); print(sort(a, decreasing=TRUE)[1:15])
cat("\nunit_of_measurement table:\n"); print(table(dat$unit_of_measurement, useNA="ifany"))
'n_after = 1338
=== altitude_mean_meters summary ===
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 1100 1311 1775 1600 190164 230
NA count: 230
zero/neg count: 0
top values:
[1] 190164 190164 110000 11000 4287 4001 3850 3845 3825 3800
[11] 3500 3280 3280 3280 3170
unit_of_measurement table:
ft m
182 1156cd /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)
a <- dat$altitude_mean_meters
pos <- a[!is.na(a) & a>0]
la <- log10(pos)
q <- quantile(la, c(.25,.75))
iqr <- q[[2]]-q[[1]]
lo <- q[[1]]-3*iqr; hi <- q[[2]]+3*iqr
cat("log10 Q1,Q3:", q[[1]], q[[2]], " IQR:", iqr, "\n")
cat("log10 fence:", lo, hi, "\n")
cat("meter fence:", 10^lo, 10^hi, "\n")
flag <- !is.na(a) & (log10(a) < lo | log10(a) > hi)
cat("n flagged (incl NA-safe):", sum(flag), "\n")
fr <- dat[flag, c("altitude","altitude_mean_meters","unit_of_measurement")]
print(as.data.frame(fr), row.names=FALSE)
'log10 Q1,Q3: 3.041393 3.20412 IQR: 0.1627273
log10 fence: 2.553211 3.692302
meter fence: 357.4463 4923.817
n flagged (incl NA-safe): 51
altitude altitude_mean_meters unit_of_measurement
1.2 12.0000 m
1.2 12.0000 m
-1 1.0000 m
350 350.0000 m
170 170.0000 m
approx. 1'500 metres a.s.l. 250.5000 m
1 1.0000 m
175 175.0000 m
16.8 168.0000 m
approx.1'500 metres a.s.l. 250.5000 m
1 1.0000 m
350 350.0000 m
11000 metros 11000.0000 m
518 157.8864 ft
1.65 165.0000 m
1 1.0000 m
200 m 200.0000 m
518 157.8864 ft
1 1.0000 m
200 m 200.0000 m
1 1.0000 m
1 1.0000 m
1 1.0000 m
1.3 13.0000 m
518 157.8864 ft
1 1.0000 m
1 1.0000 m
1 1.0000 m
150 m 150.0000 m
190164 190164.0000 m
1.3 13.0000 m
160-200m 180.0000 m
1.25 125.0000 m
110公尺 110.0000 m
1 1.0000 m
200 200.0000 m
300 m 300.0000 m
100 100.0000 m
1100.00 mosl 110000.0000 m
1.85 185.0000 m
50 50.0000 m
280 280.0000 m
200 m 200.0000 m
150 m 150.0000 m
250 m 250.0000 m
250 m 250.0000 m
1901.64 190164.0000 m
200 200.0000 m
12oo 12.0000 m
~350m 350.0000 m
40 40.0000 mcd /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)
a <- dat$altitude_mean_meters
pos <- a[!is.na(a) & a>0]
la <- log10(pos); q <- quantile(la, c(.25,.75)); iqr <- q[[2]]-q[[1]]
lo <- 10^(q[[1]]-3*iqr); hi <- 10^(q[[2]]+3*iqr)
flag <- !is.na(a) & (a < lo | a > hi)
first_tok <- function(s){ m <- str_extract(s, "[0-9]+(\\.[0-9]+)?"); as.numeric(m) }
idx <- which(flag)
corr <- rep(NA_real_, nrow(dat)); ncorr <- 0L
for(i in idx){
tok <- first_tok(dat$altitude[i])
if(is.na(tok)) next
cands <- c(tok/10, tok/100, tok)
hit <- cands[cands>=lo & cands<=hi]
if(length(hit)>0){ corr[i] <- hit[1]; ncorr <- ncorr+1L }
}
cat("fence m:", round(lo,3), round(hi,3), "\n")
cat("n flagged:", length(idx), " n corrected:", ncorr, "\n\n")
show <- dat[idx, c("altitude","altitude_mean_meters")]
show$first_tok <- sapply(dat$altitude[idx], first_tok)
show$corrected <- corr[idx]
print(as.data.frame(show), row.names=FALSE)
'fence m: 357.446 4923.817
n flagged: 51 n corrected: 7
altitude altitude_mean_meters first_tok corrected
1.2 12.0000 1.20 NA
1.2 12.0000 1.20 NA
-1 1.0000 1.00 NA
350 350.0000 350.00 NA
170 170.0000 170.00 NA
approx. 1'500 metres a.s.l. 250.5000 1.00 NA
1 1.0000 1.00 NA
175 175.0000 175.00 NA
16.8 168.0000 16.80 NA
approx.1'500 metres a.s.l. 250.5000 1.00 NA
1 1.0000 1.00 NA
350 350.0000 350.00 NA
11000 metros 11000.0000 11000.00 1100.00
518 157.8864 518.00 518.00
1.65 165.0000 1.65 NA
1 1.0000 1.00 NA
200 m 200.0000 200.00 NA
518 157.8864 518.00 518.00
1 1.0000 1.00 NA
200 m 200.0000 200.00 NA
1 1.0000 1.00 NA
1 1.0000 1.00 NA
1 1.0000 1.00 NA
1.3 13.0000 1.30 NA
518 157.8864 518.00 518.00
1 1.0000 1.00 NA
1 1.0000 1.00 NA
1 1.0000 1.00 NA
150 m 150.0000 150.00 NA
190164 190164.0000 190164.00 1901.64
1.3 13.0000 1.30 NA
160-200m 180.0000 160.00 NA
1.25 125.0000 1.25 NA
110公尺 110.0000 110.00 NA
1 1.0000 1.00 NA
200 200.0000 200.00 NA
300 m 300.0000 300.00 NA
100 100.0000 100.00 NA
1100.00 mosl 110000.0000 1100.00 1100.00
1.85 185.0000 1.85 NA
50 50.0000 50.00 NA
280 280.0000 280.00 NA
200 m 200.0000 200.00 NA
150 m 150.0000 150.00 NA
250 m 250.0000 250.00 NA
250 m 250.0000 250.00 NA
1901.64 190164.0000 1901.64 1901.64
200 200.0000 200.00 NA
12oo 12.0000 12.00 NA
~350m 350.0000 350.00 NA
40 40.0000 40.00 NAcd /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(col in c("category_one_defects","category_two_defects","quakers")){
x <- dat[[col]]
cat("=== ", col, " ===\n")
cat("NA:", sum(is.na(x)), " zeros:", sum(x==0,na.rm=TRUE), sprintf("(%.1f%%)", 100*mean(x==0,na.rm=TRUE)), "\n")
cat("max:", max(x,na.rm=TRUE), " mean:", round(mean(x,na.rm=TRUE),3), "\n")
print(quantile(x, c(.5,.75,.9,.95,.98,.99,.995,1), na.rm=TRUE))
q<-quantile(x,c(.25,.75),na.rm=TRUE); cat("raw IQR3 upper fence:", q[[2]]+3*(q[[2]]-q[[1]]), "\n\n")
}
'=== category_one_defects ===
NA: 0 zeros: 1136 (84.9%)
max: 63 mean: 0.48
50% 75% 90% 95% 98% 99% 99.5% 100%
0.000 0.000 1.000 2.000 5.000 7.630 11.315 63.000
raw IQR3 upper fence: 0
=== category_two_defects ===
NA: 0 zeros: 373 (27.9%)
max: 55 mean: 3.558
50% 75% 90% 95% 98% 99% 99.5% 100%
2.000 4.000 9.000 13.000 20.000 26.630 31.315 55.000
raw IQR3 upper fence: 16
=== quakers ===
NA: 1 zeros: 1243 (93.0%)
max: 11 mean: 0.174
50% 75% 90% 95% 98% 99% 99.5% 100%
0 0 0 1 2 5 6 11
raw IQR3 upper fence: 0cd /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(col in c("category_one_defects","category_two_defects","quakers")){
x <- dat[[col]]; xnn <- x[!is.na(x)]
pos <- xnn[xnn>0]
# Tukey 1.5 on positive-only
qp<-quantile(pos,c(.25,.75)); up_pos<-qp[[2]]+1.5*(qp[[2]]-qp[[1]])
# percentile thresholds on full col
q98<-quantile(xnn,.98); q99<-quantile(xnn,.99)
cat(sprintf("%-22s pos-Tukey1.5 up=%.2f n>%.2f=%d (%.2f%%) | q98=%.1f n=%d (%.2f%%) | q99=%.1f n=%d (%.2f%%)\n",
col, up_pos, up_pos, sum(xnn>up_pos), 100*mean(xnn>up_pos),
q98, sum(xnn>q98), 100*mean(xnn>q98),
q99, sum(xnn>q99), 100*mean(xnn>q99)))
}
'category_one_defects pos-Tukey1.5 up=6.00 n>6.00=17 (1.27%) | q98=5.0 n=20 (1.49%) | q99=7.6 n=14 (1.05%) category_two_defects pos-Tukey1.5 up=12.00 n>12.00=73 (5.46%) | q98=20.0 n=23 (1.72%) | q99=26.6 n=14 (1.05%) quakers pos-Tukey1.5 up=6.00 n>6.00=6 (0.45%) | q98=2.0 n=25 (1.87%) | q99=5.0 n=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)
prim <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
clean <- c("uniformity","clean_cup","sweetness")
cat("PRIMARY 7:\n")
for(c in prim) cat(sprintf(" %-14s mean=%.2f sd=%.3f min=%.2f\n",c,mean(dat[[c]]),sd(dat[[c]]),min(dat[[c]])))
cat("CLEANLINESS (near-constant):\n")
for(c in clean) cat(sprintf(" %-14s mean=%.2f sd=%.3f min=%.2f frac@10=%.2f\n",c,mean(dat[[c]]),sd(dat[[c]]),min(dat[[c]]),mean(dat[[c]]==10)))
G <- as.matrix(dat[,prim]); cc<-complete.cases(G)
cat("\ncomplete cases on 7 primary:", sum(cc), "of", nrow(dat),"\n")
cat("chisq 0.999 df=7 threshold:", qchisq(0.999,7), "\n")
md2<-mahalanobis(G[cc,],colMeans(G[cc,]),cov(G[cc,]))
cat("n flagged (md2>thr):", sum(md2>qchisq(0.999,7)), "\n")
cat("condition number cov(7):", kappa(cov(G[cc,])), "\n")
cat("condition number cov(10):", kappa(cov(as.matrix(dat[complete.cases(dat[,c(prim,clean)]),c(prim,clean)]))), "\n")
'PRIMARY 7: aroma mean=7.57 sd=0.316 min=5.08 flavor mean=7.53 sd=0.341 min=6.08 aftertaste mean=7.41 sd=0.350 min=6.17 acidity mean=7.54 sd=0.319 min=5.25 body mean=7.52 sd=0.308 min=5.08 balance mean=7.52 sd=0.354 min=5.25 cupper_points mean=7.51 sd=0.427 min=5.17 CLEANLINESS (near-constant): uniformity mean=9.84 sd=0.485 min=6.00 frac@10=0.86 clean_cup mean=9.84 sd=0.715 min=0.00 frac@10=0.91 sweetness mean=9.86 sd=0.554 min=1.33 frac@10=0.91 complete cases on 7 primary: 1338 of 1338 chisq 0.999 df=7 threshold: 24.32189 n flagged (md2>thr): 40 condition number cov(7): 22.05133 condition number cov(10): 69.38935
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)
prim <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
# altitude corrected vector (from earlier logic)
a <- dat$altitude_mean_meters; pos<-a[!is.na(a)&a>0]; la<-log10(pos)
q<-quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]; lo<-10^(q[[1]]-3*iqr); hi<-10^(q[[2]]+3*iqr)
flag<-!is.na(a)&(a<lo|a>hi)
ft<-function(s){as.numeric(str_extract(s,"[0-9]+(\\.[0-9]+)?"))}
alt_eff <- a
for(i in which(flag)){tok<-ft(dat$altitude[i]); if(is.na(tok))next; cd<-c(tok/10,tok/100,tok); h<-cd[cd>=lo&cd<=hi]; if(length(h)>0) alt_eff[i]<-h[1]}
cc <- complete.cases(dat[,prim]) # all TRUE
fr <- dat[cc,]
alt_eff_fr <- alt_eff[cc]
# model frames
log_alt <- log10(alt_eff_fr); log_alt[!is.finite(log_alt)] <- NA
log_alt[is.na(log_alt)] <- mean(log_alt, na.rm=TRUE)
d1<-fr$category_one_defects; d2<-fr$category_two_defects; dq<-fr$quakers
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
l1<-log1p(imp(d1)); l2<-log1p(imp(d2)); lq<-log1p(imp(dq))
base <- fr[,c("total_cup_points",prim)]
m1 <- lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points, data=base)
m2 <- lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+log_alt, data=cbind(base,log_alt))
m3 <- lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1+l2+lq, data=cbind(base,l1,l2,lq))
cat("AIC grades_only :", AIC(m1), "\n")
cat("AIC grades_altitude :", AIC(m2), "\n")
cat("AIC grades_defects :", AIC(m3), "\n")
cat("n obs:", nobs(m1), nobs(m2), nobs(m3), "\n")
cat("R2:", summary(m1)$r.squared, summary(m2)$r.squared, summary(m3)$r.squared, "\n")
'AIC grades_only : 4459.362 AIC grades_altitude : 4461.292 AIC grades_defects : 4418.68 n obs: 1338 1338 1338 R2: 0.7756335 0.7756451 0.7833264
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)
prim <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
l1<-log1p(imp(dat$category_one_defects)); l2<-log1p(imp(dat$category_two_defects)); lq<-log1p(imp(dat$quakers))
base <- cbind(dat[,c("total_cup_points",prim)],l1,l2,lq)
m3 <- lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1+l2+lq, data=base)
n<-nobs(m3); p<-length(coef(m3))
hat<-hatvalues(m3); cook<-cooks.distance(m3)
lev_t<-2*p/n; cook_t<-4/n
cat("n=",n," p=",p," lev_t=",lev_t," cook_t=",cook_t,"\n")
cat("n_high_lev:",sum(hat>lev_t)," n_high_cook:",sum(cook>cook_t)," n_both:",sum(hat>lev_t & cook>cook_t),"\n")
top<-order(-cook)[1:10]
cat("top10 cook row_ids (0-indexed):", (top-1), "\n")
cat("top10 cook values:", round(cook[top],4),"\n")
'n= 1338 p= 11 lev_t= 0.01644245 cook_t= 0.002989537 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)})
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) %>%
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_countries=%d n_changes>=2=%d max_chg=%d\n", tr, nrow(ct), sum(ct$chg>=2), max(ct$chg)))
}
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=0.1), .groups="drop") %>%
filter(n>=10) %>% mutate(raw_rank=rank(-raw_mean,ties.method="min"), trimmed_rank=rank(-trimmed_mean,ties.method="min"), chg=abs(raw_rank-trimmed_rank)) %>% arrange(raw_rank)
print(as.data.frame(ct), row.names=FALSE, digits=5)
'trim=0.1: n_countries=21 n_changes>=2=2 max_chg=2
trim=0.2: n_countries=21 n_changes>=2=2 max_chg=2
country n raw_mean trimmed_mean raw_rank trimmed_rank
Ethiopia 44 85.484 85.518 1 1
United States 10 84.433 84.905 2 2
Kenya 25 84.310 84.508 3 3
Uganda 36 83.452 83.448 4 4
Colombia 183 83.107 83.247 5 5
El Salvador 21 83.053 83.109 6 6
China 16 82.927 82.971 7 8
Costa Rica 51 82.789 83.020 8 7
Thailand 32 82.574 82.619 9 10
Indonesia 20 82.566 82.774 10 9
Peru 10 82.526 82.439 11 12
Brazil 132 82.406 82.513 12 11
Tanzania, United Republic Of 40 82.370 82.272 13 13
Taiwan 75 82.001 81.998 14 16
Guatemala 181 81.847 82.168 15 14
United States (Hawaii) 73 81.820 82.076 16 15
Malawi 11 81.712 81.712 17 17
India 14 81.083 81.382 18 19
Mexico 236 80.890 81.246 19 20
Honduras 52 80.884 81.542 20 18
Nicaragua 26 80.458 80.894 21 21
chg
0
0
0
0
0
0
1
1
1
1
1
1
0
2
1
1
0
1
1
2
0cd /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)
prim <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
# altitude flag
a<-dat$altitude_mean_meters; pos<-a[!is.na(a)&a>0]; la<-log10(pos)
q<-quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]; lo<-10^(q[[1]]-3*iqr); hi<-10^(q[[2]]+3*iqr)
flag_alt<-!is.na(a)&(a<lo|a>hi)
# mv flag
G<-as.matrix(dat[,prim]); md2<-mahalanobis(G,colMeans(G),cov(G)); flag_mv<-md2>qchisq(0.999,7)
# influence on grades_defects
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
l1<-log1p(imp(dat$category_one_defects)); l2<-log1p(imp(dat$category_two_defects)); lq<-log1p(imp(dat$quakers))
base<-cbind(dat[,c("total_cup_points",prim)],l1,l2,lq)
m3<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1+l2+lq,data=base)
n<-nobs(m3);p<-length(coef(m3)); hat<-hatvalues(m3);cook<-cooks.distance(m3)
hl<-hat>2*p/n; hc<-cook>4/n
composite <- flag_alt | flag_mv | (hl & hc)
cat("n composite:", sum(composite), "\n")
# sensitivity point est
f_full<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat)
f_drop<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat[!composite,])
cf<-unname(coef(f_full)["flavor"]); cd<-unname(coef(f_drop)["flavor"])
cat("coef full:",cf," dropped:",cd," delta_pct:",100*(cd-cf)/cf,"\n")
# bootstrap
set.seed(20260512); B<-500; N<-nrow(dat); deltas<-numeric(0)
for(b in 1:B){
idx<-sample.int(N,N,replace=TRUE)
db<-dat[idx,]; cb<-composite[idx]
kept<-db[!cb,]; if(nrow(kept)<10) next
ff<-tryCatch(lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=db),error=function(e)NULL)
fd<-tryCatch(lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=kept),error=function(e)NULL)
if(is.null(ff)||is.null(fd))next
cff<-unname(coef(ff)["flavor"]); cfd<-unname(coef(fd)["flavor"])
if(is.na(cff)||is.na(cfd)||cff==0)next
deltas<-c(deltas,100*(cfd-cff)/cff)
}
ci<-quantile(deltas,c(.025,.975),names=FALSE)
cat("n_boot:",length(deltas)," CI:",ci[1],ci[2],"\n")
'n composite: 105 coef full: 2.085398 dropped: 1.652728 delta_pct: -20.74764 n_boot: 500 CI: -44.32811 6.247039
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)
prim <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
a<-dat$altitude_mean_meters; pos<-a[!is.na(a)&a>0]; la<-log10(pos)
q<-quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]; lo<-10^(q[[1]]-3*iqr); hi<-10^(q[[2]]+3*iqr)
A<-!is.na(a)&(a<lo|a>hi)
G<-as.matrix(dat[,prim]); M<-mahalanobis(G,colMeans(G),cov(G))>qchisq(0.999,7)
q98<-function(x)quantile(x,.98,na.rm=TRUE)
D<-(!is.na(dat$category_one_defects)&dat$category_one_defects>q98(dat$category_one_defects)) |
(!is.na(dat$category_two_defects)&dat$category_two_defects>q98(dat$category_two_defects)) |
(!is.na(dat$quakers)&dat$quakers>q98(dat$quakers))
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
l1<-log1p(imp(dat$category_one_defects)); l2<-log1p(imp(dat$category_two_defects)); lq<-log1p(imp(dat$quakers))
base<-cbind(dat[,c("total_cup_points",prim)],l1,l2,lq)
m3<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1+l2+lq,data=base)
n<-nobs(m3);p<-length(coef(m3)); I<-(hatvalues(m3)>2*p/n)&(cooks.distance(m3)>4/n)
tot<-as.integer(A)+as.integer(M)+as.integer(D)+as.integer(I)
cat("flags: A",sum(A)," M",sum(M)," D",sum(D)," I",sum(I),"\n")
cat("0:",sum(tot==0)," 1:",sum(tot==1)," 2:",sum(tot==2)," 3:",sum(tot==3)," 4:",sum(tot==4),"\n")
cat("alt_only:",sum(A&!M&!I&!D)," mv_only:",sum(M&!A&!I&!D)," infl_only:",sum(I&!A&!M&!D)," def_only:",sum(D&!A&!M&!I),"\n")
cat("alt&mv:",sum(A&M)," alt&infl:",sum(A&I)," mv&infl:",sum(M&I)," alt&mv&infl:",sum(A&M&I),"\n")
cat("sum check:",sum(tot==0)+sum(tot==1)+sum(tot==2)+sum(tot==3)+sum(tot==4),"=",nrow(dat),"\n")
'flags: A 51 M 40 D 63 I 35 0: 1186 1: 122 2: 23 3: 7 4: 0 alt_only: 42 mv_only: 23 infl_only: 10 def_only: 47 alt&mv: 4 alt&infl: 5 mv&infl: 13 alt&mv&infl: 1 sum check: 1338 = 1338
command: create
path: /app/analysis.R
file_text: #!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-matched rebuild
# =====================================================
# Each column is audited with a method matched to its distribution rather than
# one off-the-shelf recipe: log-scale fences for right-skewed altitude (plus
# unit-slip repair), tail-quantile thresholds for zero-inflated defect counts,
# classical Mahalanobis on the seven informative grades only, AIC-selected OLS
# influence (leverage AND Cook), and a robust (trimmed-mean) country ranking
# with a composite-drop sensitivity bootstrap.
suppressPackageStartupMessages({
library(dplyr)
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)
PRIM <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "cupper_points")
# ----------------------------------------------------------------------------
# 1. Load + sentinel drop
# ----------------------------------------------------------------------------
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
# Withdrawn submission: every grade recorded as zero. Drop before anything else.
zero_score <- with(raw, aroma == 0 & flavor == 0 & aftertaste == 0 &
acidity == 0 & body == 0 & balance == 0 & cupper_points == 0)
zero_score[is.na(zero_score)] <- FALSE
dat <- raw[!zero_score, , drop = FALSE]
n_after <- nrow(dat)
dat$row_id <- as.integer(seq_len(n_after) - 1L)
# ----------------------------------------------------------------------------
# 2. Altitude: log10 Tukey fence (k = 3) + decimal-displacement repair
# ----------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
la <- log10(alt_pos)
qa <- quantile(la, c(0.25, 0.75), names = FALSE)
iqr_la <- qa[2] - qa[1]
lo_log <- qa[1] - 3 * iqr_la
hi_log <- qa[2] + 3 * iqr_la
lo_m <- 10^lo_log
hi_m <- 10^hi_log
flag_alt <- !is.na(alt) & (alt < lo_m | alt > hi_m)
n_flag_alt <- sum(flag_alt)
# Try to recover unit slips on the raw altitude string: first numeric token,
# tested as /10, then /100, then as-is; keep the first candidate inside fence.
first_token <- function(s) as.numeric(str_extract(s, "[0-9]+(\\.[0-9]+)?"))
alt_corrected <- rep(NA_real_, n_after) # reported only when a slip is fixed
alt_effective <- alt # value fed to the regression
n_unit_corrected <- 0L
for (i in which(flag_alt)) {
tok <- first_token(dat$altitude[i])
if (is.na(tok)) next
cands <- c(tok / 10, tok / 100, tok)
inb <- cands[cands >= lo_m & cands <= hi_m]
if (length(inb) > 0) {
alt_corrected[i] <- inb[1]
alt_effective[i] <- inb[1]
n_unit_corrected <- n_unit_corrected + 1L
}
}
# ----------------------------------------------------------------------------
# 3. Defects: tail-quantile thresholds (zero-inflated counts)
# ----------------------------------------------------------------------------
# Raw IQR fences collapse to 0 when Q1 = Q3 = 0, so flag the sparse upper tail
# via the 98th empirical percentile (strictly above) -> ~top 1-2% per column.
defect_tail <- function(x) {
thr <- unname(quantile(x, 0.98, na.rm = TRUE, type = 7))
fl <- !is.na(x) & x > thr
list(upper_threshold = thr, n = sum(fl), flag = fl)
}
d1 <- defect_tail(dat$category_one_defects)
d2 <- defect_tail(dat$category_two_defects)
dq <- defect_tail(dat$quakers)
flag_defect <- d1$flag | d2$flag | dq$flag
# ----------------------------------------------------------------------------
# 4. Multivariate grades: classical Mahalanobis on the 7 informative grades
# ----------------------------------------------------------------------------
G <- as.matrix(dat[, PRIM])
cc_g <- complete.cases(G)
mu <- colMeans(G[cc_g, , drop = FALSE])
S <- cov(G[cc_g, , drop = FALSE])
md2
…[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 altitude fence [357.45, 4923.82] m; flagged=51 corrected=7 mv flagged=40; influence both=35; composite=105 flavor coef 2.0854 -> 1.6527 (-20.75%); boot CI [-44.33, 6.25] n=500
cd /app && cat outputs/outlier_report.json
{
"n_input": 1339,
"n_after_zero_score_drop": 1338,
"methods": {
"altitude": "log10 Tukey fence (k=3) on positive altitude_mean_meters, back-transformed to meters; flagged rows repaired by testing first-token /10, /100, as-is against the fence",
"defects": "per-column 98th-percentile upper-tail threshold on zero-inflated counts (flag x>thr), since raw IQR fences collapse to 0",
"multivariate_grades": "classical Mahalanobis (sample mean/cov) on the 7 informative SCA grades; chi-square 99.9% df=7 cutoff",
"influence": "AIC-selected OLS among grades_only/grades_altitude/grades_defects; flag leverage>2p/n AND Cook's D>4/n on the selected fit",
"country_ranking": "countries with n>=10: raw mean vs 10% trimmed mean, compared by descending rank (ties=min)"
},
"altitude": {
"lower_bound_m": 357.446289062499,
"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.3218863478569,
"n_flagged": 40
},
"influence": {
"selected_model": "grades_defects",
"aic_grades_only": 4459.36151596946,
"aic_grades_altitude": 4461.29232958899,
"aic_grades_defects": 4418.68007429542,
"leverage_threshold": 0.0164424514200299,
"cook_threshold": 0.00298953662182362,
"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.4840909090909,
"trimmed_mean": 85.5177777777778,
"raw_rank": 1,
"trimmed_rank": 1
},
{
"country": "United States",
"n": 10,
"raw_mean": 84.433,
"trimmed_mean": 84.905,
"raw_rank": 2,
"trimmed_rank": 2
},
{
"country": "Kenya",
"n": 25,
"raw_mean": 84.3096,
"trimmed_mean": 84.507619047619,
"raw_rank": 3,
"trimmed_rank": 3
},
{
"country": "Uganda",
"n": 36,
"raw_mean": 83.4519444444445,
"trimmed_mean": 83.448,
"raw_rank": 4,
"trimmed_rank": 4
},
{
"country": "Colombia",
"n": 183,
"raw_mean": 83.1065573770492,
"trimmed_mean": 83.2466666666667,
"raw_rank": 5,
"trimmed_rank": 5
},
{
"country": "El Salvador",
"n": 21,
"raw_mean": 83.0528571428571,
"trimmed_mean": 83.1094117647059,
"raw_rank": 6,
"trimmed_rank": 6
},
{
"country": "China",
"n": 16,
"raw_mean": 82.9275,
"trimmed_mean": 82.9707142857143,
"raw_rank": 7,
"trimmed_rank": 8
},
{
"country": "Costa Rica",
"n": 51,
"raw_mean": 82.7890196078431,
"trimmed_mean": 83.02,
"raw_rank": 8,
"trimmed_rank": 7
},
{
"country": "Thailand",
"n": 32,
"raw_mean": 82.57375,
"trimmed_mean": 82.6192307692308,
"raw_rank": 9,
"trimmed_rank": 10
},
{
"country": "Indonesia",
"n": 20,
"raw_mean": 82.5655,
"trimmed_mean": 82.77375,
"raw_rank": 10,
"trimmed_rank": 9
},
{
"country": "Peru",
"n": 10,
"raw_mean": 82.526,
"trimmed_mean": 82.43875,
"raw_rank": 11,
"trimmed_rank": 12
}
…[truncated]cd /app && python3 -c "import json; d=json.load(open('outputs/outlier_report.json')); print('JSON valid. top keys:', list(d.keys()))" && echo "---" && Rscript -e '
suppressPackageStartupMessages({library(readr)})
fit <- readRDS("outputs/coffee_ols_selected.rds")
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
cat("class(fit):", class(fit), " formula:", deparse(formula(fit)), "\n")
n<-nobs(fit); p<-length(coef(fit)); cat("n=",n," p=",p,"\n")
hat<-hatvalues(fit); cook<-cooks.distance(fit)
pos<-as.integer(rownames(model.frame(fit)))
hl<-rep(FALSE,nrow(csv)); hc<-rep(FALSE,nrow(csv))
hl[pos]<-hat>2*p/n; hc[pos]<-cook>4/n
cat("verifier hl==csv:", all(hl==csv$high_leverage_flag), " hc==csv:", all(hc==csv$high_cook_flag),"\n")
cat("n rows csv:", nrow(csv), " row_id range:", min(csv$row_id), max(csv$row_id),"\n")
comp_check <- csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag | (csv$high_leverage_flag & csv$high_cook_flag)
cat("composite reconciles:", all(comp_check==csv$composite_drop_flag),"\n")
# top10 cook check
co<-rep(NA_real_,nrow(csv)); co[pos]<-cook
cat("top10 ids:", paste((order(co,decreasing=TRUE,na.last=NA)[1:10])-1,collapse=" "),"\n")
cat("corrected non-NA count:", sum(!is.na(csv$altitude_corrected_m)), "\n")
cat("corrected rows all altitude-flagged:", all(csv$altitude_outlier_flag[!is.na(csv$altitude_corrected_m)]),"\n")
'JSON valid. top keys: ['n_input', 'n_after_zero_score_drop', 'methods', 'altitude', 'defects', 'multivariate_grades', 'influence', 'country_ranking', 'sensitivity', 'flag_interactions'] --- class(fit): lm formula: total_cup_points ~ aroma + flavor + aftertaste + acidity + body + balance + cupper_points + l1 + l2 + lq n= 1338 p= 11 verifier hl==csv: TRUE hc==csv: TRUE n rows csv: 1338 row_id range: 0 1337 composite reconciles: TRUE top10 ids: 1309 1308 1333 1295 1306 1260 1307 1302 1291 1259 corrected non-NA count: 7 corrected rows all altitude-flagged: TRUE
cd /app && head -3 outputs/outlier_flags.csv && echo "---CORRECTED ROWS---" && Rscript -e '
suppressPackageStartupMessages(library(readr))
csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
raw<-read_csv("data/coffee_ratings.csv",show_col_types=FALSE)
dat<-raw[!(raw$aroma==0&raw$flavor==0&raw$aftertaste==0&raw$acidity==0&raw$body==0&raw$balance==0&raw$cupper_points==0),]
sub<-csv[!is.na(csv$altitude_corrected_m),]
sub$raw_altitude<-dat$altitude[sub$row_id+1]
sub$alt_mean<-dat$altitude_mean_meters[sub$row_id+1]
print(as.data.frame(sub[,c("row_id","raw_altitude","alt_mean","altitude_corrected_m")]),row.names=FALSE)
cat("\ncol names:",paste(names(csv),collapse=", "),"\n")
cat("logical cols are TRUE/FALSE text:\n"); print(sapply(csv[,c("altitude_outlier_flag","defect_outlier_flag","composite_drop_flag")],class))
'row_id,country_of_origin,total_cup_points,altitude_outlier_flag,altitude_corrected_m,defect_outlier_flag,multivariate_grade_outlier_flag,high_leverage_flag,high_cook_flag,composite_drop_flag
0,Ethiopia,90.58,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE
1,Ethiopia,89.92,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE
---CORRECTED ROWS---
row_id raw_altitude alt_mean altitude_corrected_m
543 11000 metros 11000.0000 1100.00
584 518 157.8864 518.00
677 518 157.8864 518.00
833 518 157.8864 518.00
896 190164 190164.0000 1901.64
1040 1100.00 mosl 110000.0000 1100.00
1144 1901.64 190164.0000 1901.64
col names: row_id, country_of_origin, total_cup_points, altitude_outlier_flag, altitude_corrected_m, defect_outlier_flag, multivariate_grade_outlier_flag, high_leverage_flag, high_cook_flag, composite_drop_flag
logical cols are TRUE/FALSE text:
altitude_outlier_flag defect_outlier_flag composite_drop_flag
"logical" "logical" "logical"command: create path: /app/outputs/audit_memo.md file_text: # Coffee Quality Outlier Audit , Method Memo The draft applied one off-the-shelf recipe (Tukey fences, Mahalanobis on every grade, `|rstandard|>2`, raw country means) to columns with very different shapes. Each block below states why that default is wrong for the column and what shape-matched method replaced it. ## Data prep One lot (Honduras) is a withdrawn submission with **every grade recorded as zero**; its `total_cup_points` is 0. The draft's `total_cup_points > 0` filter happens to remove it, but only as a side effect. I drop it explicitly on the all-grades-zero sentinel and report counts: **n_input = 1339**, **n_after_zero_score_drop = 1338**. All 1338 survivors are complete on the seven primary grades, so downstream frames align one-to-one with `row_id`. ## Altitude `altitude_mean_meters` is strongly right-skewed (median 1311 m, max 190,164 m). Tukey fences on the **raw meter scale** are driven by that skew and are meaningless here. I take `log10` over positive values, build the IQR fence with `k = 3` on the log scale, and back-transform: the fence is **[357 m, 4924 m]**, flagging **51** rows. Many flags are decimal-displacement typos in the raw `altitude` string (`190164`, `11000 metros`, `1100.00 mosl`). For each flagged row I take the first numeric token and test `/10`, `/100`, then as-is, keeping the first candidate inside the fence; this repairs **7** rows (e.g. `190164 → 1901.64`, `11000 → 1100`). Rows with no in-bounds candidate (genuine lowland lots coded `1`, `50`, `200 m`) keep `altitude_corrected_m = NA` rather than being silently rescaled. The draft never attempted any repair. ## Defects `category_one_defects`, `category_two_defects`, and `quakers` are mass-at-zero (85%, 28%, 93% zeros). For cat-one and quakers `Q1 = Q3 = 0`, so the draft's raw IQR fence collapses to `upper = 0` and would flag **every** non-zero lot , hundreds of rows, not outliers. I instead threshold each column at its **98th empirical percentile** and flag counts strictly above it, which isolates the sparse upper tail the prompt describes: thresholds 5 / 20 / 2 flagging 20 / 23 / 25 rows (1.5% / 1.7% / 1.9%). A lot is a defect outlier if any column trips. ## Multivariate grades The draft ran Mahalanobis on **all ten** grade columns, including `uniformity`, `clean_cup`, and `sweetness`, which are near-constant (≈90% pinned at 10). Those spikes inflate the covariance condition number (≈69 vs ≈22) and make the chi-squared reference distribution invalid. I restrict to the **seven informative SCA attributes** (`aroma`, `flavor`, `aftertaste`, `acidity`, `body`, `balance`, `cupper_points`) and use a classical Mahalanobis distance (sample mean and covariance) against a **χ²(0.999, df = 7) = 24.32** cutoff, flagging **40** joint outliers. MCD is not required here; the classical estimator with an explicit high quantile is what the spec asks for. ## Influence `|rstandard| > 2` measures residual size, not **influence**, and flags ~5% of any clean model by construction; the draft also reused that one vector for leverage, Cook, and "both". I fit three OLS candidates on the 1338 complete-grade rows and pick the lowest finite AIC: grades_only (4459.4), grades_altitude (4461.3, adds mean-imputed `log10` elevation using corrected meters), and **grades_defects (4418.7, selected)**, which adds `log1p` of the three defect columns with mean-imputed NAs. All three share the same response and n, so AIC is comparable. On the selected fit (n = 1338, p = 11) I flag **high leverage** (`hat > 2p/n = 0.0164`, 110 rows) and **high Cook's D** (`> 4/n = 0.00299`, 66 rows) separately; the influence flag is their **intersection** (35 rows). The fitted model is saved to `coffee_ols_selected.rds` so the verifier can recompute both diagnostics. ## Country ranking Raw country means let a single bad lot swing a country's standing and offer no robustness check. For the **21** countries with ≥10 lots I compare t …[truncated]
New file created successfully at /app/outputs/audit_memo.md
cd /app && echo "Word count: $(wc -w < outputs/audit_memo.md)" && echo "---Headings---" && grep -E '^## ' outputs/audit_memo.md && echo "---Required check---" && for h in "Data prep" "Altitude" "Defects" "Multivariate grades" "Influence" "Country ranking" "Sensitivity" "Composite flags"; do grep -q "^## $h\$" outputs/audit_memo.md && echo "OK: $h" || echo "MISSING: $h"; done
Word count: 844 ---Headings--- ## Data prep ## Altitude ## Defects ## Multivariate grades ## Influence ## Country ranking ## Sensitivity ## Composite flags ---Required check--- OK: Data prep OK: Altitude OK: Defects OK: Multivariate grades OK: Influence OK: Country ranking OK: Sensitivity OK: Composite flags
command: str_replace
path: /app/analysis.R
old_str: # ----------------------------------------------------------------------------
# Write per-row flags CSV
new_str: # ----------------------------------------------------------------------------
# Descriptive stats used only for the memo prose (kept reproducible)
# ----------------------------------------------------------------------------
z_frac <- function(x) 100 * mean(x == 0, na.rm = TRUE)
clean_at10 <- mean(vapply(c("uniformity", "clean_cup", "sweetness"),
function(c) mean(dat[[c]] == 10, na.rm = TRUE), numeric(1)))
cond_7 <- kappa(S)
cond_10 <- {
g10 <- c(PRIM, "uniformity", "clean_cup", "sweetness")
M10 <- as.matrix(dat[complete.cases(dat[, g10]), g10])
kappa(cov(M10))
}
movers <- country_tbl %>% filter(abs_rank_change >= 2)
movers_str <- if (nrow(movers) > 0)
paste(sprintf("%s %d\u2192%d", movers$country, movers$raw_rank, movers$trimmed_rank),
collapse = ", ") else "none"
pct <- function(k) 100 * k / n_after
# ----------------------------------------------------------------------------
# Write per-row flags CSVThe file /app/analysis.R has been edited successfully.
command: str_replace
path: /app/analysis.R
old_str: paste(sprintf("%s %d\u2192%d", movers$country, movers$raw_rank, movers$trimmed_rank),
collapse = ", ") else "none"
new_str: paste(sprintf("%s %d->%d", movers$country, movers$raw_rank, movers$trimmed_rank),
collapse = ", ") else "none"The file /app/analysis.R has been edited successfully.
command: str_replace
path: /app/analysis.R
old_str: saveRDS(sel_fit, file.path(OUTPUT_DIR, "coffee_ols_selected.rds"))
new_str: saveRDS(sel_fit, file.path(OUTPUT_DIR, "coffee_ols_selected.rds"))
# ----------------------------------------------------------------------------
# Write audit_memo.md (generated here so analysis.R reproduces every artifact)
# ----------------------------------------------------------------------------
memo <- c(
"# Coffee Quality Outlier Audit - Method Memo",
"",
"The draft applied one off-the-shelf recipe (Tukey fences, Mahalanobis on every",
"grade, `|rstandard|>2`, raw country means) to columns with very different",
"shapes. Each block states why that default is wrong for the column and what",
"shape-matched method replaced it.",
"",
"## Data prep",
sprintf("One lot (Honduras) is a withdrawn submission with **every grade recorded as zero** and `total_cup_points` 0. The draft's `total_cup_points > 0` filter removes it only as a side effect; I drop it explicitly on the all-grades-zero sentinel and report counts: **n_input = %d**, **n_after_zero_score_drop = %d**. All survivors are complete on the seven primary grades, so frames align one-to-one with `row_id`.", n_input, n_after),
"",
"## Altitude",
sprintf("`altitude_mean_meters` is strongly right-skewed (median ~1311 m, max 190,164 m), so raw-meter Tukey fences are driven by the skew and meaningless. I take `log10` of positive values, build the IQR fence with `k = 3` on the log scale, and back-transform: the fence is **[%.0f m, %.0f m]**, flagging **%d** rows. Many flags are decimal-displacement typos in the raw `altitude` string (`190164`, `11000 metros`, `1100.00 mosl`). For each flagged row I take the first numeric token and test `/10`, `/100`, then as-is, keeping the first candidate inside the fence; this repairs **%d** rows (e.g. 190164 -> 1901.64, 11000 -> 1100). Rows with no in-bounds candidate (genuine lowland lots coded 1, 50, 200 m) keep `altitude_corrected_m = NA` instead of being silently rescaled. The draft attempted no repair.", lo_m, hi_m, n_flag_alt, n_unit_corrected),
"",
"## Defects",
sprintf("`category_one_defects`, `category_two_defects`, and `quakers` are mass-at-zero (%.0f%%, %.0f%%, %.0f%% zeros). For cat-one and quakers Q1 = Q3 = 0, so the draft's raw IQR fence collapses to `upper = 0` and would flag **every** non-zero lot (hundreds of rows). I instead threshold each column at its **98th empirical percentile** and flag counts strictly above it, isolating the sparse upper tail: thresholds %g / %g / %g flagging %d / %d / %d rows (%.1f%% / %.1f%% / %.1f%%). A lot is a defect outlier if any column trips.", z_frac(dat$category_one_defects), z_frac(dat$category_two_defects), z_frac(dat$quakers), d1$upper_threshold, d2$upper_threshold, dq$upper_threshold, d1$n, d2$n, dq$n, pct(d1$n), pct(d2$n), pct(dq$n)),
"",
"## Multivariate grades",
sprintf("The draft ran Mahalanobis on **all ten** grade columns, including `uniformity`, `clean_cup`, and `sweetness`, which are near-constant (~%.0f%% pinned at 10). Those spikes inflate the covariance condition number (%.0f vs %.0f) and invalidate the chi-squared reference. I restrict to the **seven informative SCA attributes** (aroma, flavor, aftertaste, acidity, body, balance, cupper_points) and use a classical Mahalanobis distance (sample mean/cov) against a **chi-square(0.999, df = 7) = %.2f** cutoff, flagging **%d** joint outliers. A high-breakdown MCD is not required; the spec asks for the classical estimator with an explicit upper-tail quantile.", 100 * clean_at10, cond_10, cond_7, mv_thresh, n_flag_mv),
"",
"## Influence",
sprintf("`|rstandard| > 2` measures residual size, not **influence**, and flags ~5%% of any clean model by construction; the draft also reused that one vector for leverage, Cook, and \"both\". I fit three OLS candidates on the complete-grade rows and pick the lowest finite AIC: grades_only (%.1f), grades_altitude (%.1f, mean-imputed log10 elevation using corrected meters
…[truncated]The file /app/analysis.R has been edited successfully.
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R && echo "=== FILES ===" && ls -la outputs/ && echo "=== MEMO WORDS: $(wc -w < outputs/audit_memo.md) ===" && echo "=== HEADINGS ===" && grep -E '^## ' outputs/audit_memo.md
Audit complete: n_input=1339 n_after=1338 selected=grades_defects altitude fence [357.45, 4923.82] m; flagged=51 corrected=7 mv flagged=40; influence both=35; composite=105 flavor coef 2.0854 -> 1.6527 (-20.75%); boot CI [-44.33, 6.25] n=500 === FILES === total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:24 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:23 .. -rw-r--r-- 1 h2tagent h2tagent 5143 Jun 20 18:24 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137130 Jun 20 18:24 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78269 Jun 20 18:24 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6634 Jun 20 18:24 outlier_report.json === MEMO WORDS: 796 === === HEADINGS === ## Data prep ## Altitude ## Defects ## Multivariate grades ## Influence ## Country ranking ## Sensitivity ## Composite flags
cd /app && cat outputs/audit_memo.md
# Coffee Quality Outlier Audit - Method Memo The draft applied one off-the-shelf recipe (Tukey fences, Mahalanobis on every grade, `|rstandard|>2`, raw country means) to columns with very different shapes. Each block states why that default is wrong for the column and what shape-matched method replaced it. ## Data prep One lot (Honduras) is a withdrawn submission with **every grade recorded as zero** and `total_cup_points` 0. The draft's `total_cup_points > 0` filter removes it only as a side effect; I drop it explicitly on the all-grades-zero sentinel and report counts: **n_input = 1339**, **n_after_zero_score_drop = 1338**. All survivors are complete on the seven primary grades, so frames align one-to-one with `row_id`. ## Altitude `altitude_mean_meters` is strongly right-skewed (median ~1311 m, max 190,164 m), so raw-meter Tukey fences are driven by the skew and meaningless. I take `log10` of positive values, build the IQR fence with `k = 3` on the log scale, and back-transform: the fence is **[357 m, 4924 m]**, flagging **51** rows. Many flags are decimal-displacement typos in the raw `altitude` string (`190164`, `11000 metros`, `1100.00 mosl`). For each flagged row I take the first numeric token and test `/10`, `/100`, then as-is, keeping the first candidate inside the fence; this repairs **7** rows (e.g. 190164 -> 1901.64, 11000 -> 1100). Rows with no in-bounds candidate (genuine lowland lots coded 1, 50, 200 m) keep `altitude_corrected_m = NA` instead of being silently rescaled. The draft attempted no repair. ## Defects `category_one_defects`, `category_two_defects`, and `quakers` are mass-at-zero (85%, 28%, 93% zeros). For cat-one and quakers Q1 = Q3 = 0, so the draft's raw IQR fence collapses to `upper = 0` and would flag **every** non-zero lot (hundreds of rows). I instead threshold each column at its **98th empirical percentile** and flag counts strictly above it, isolating the sparse upper tail: thresholds 5 / 20 / 2 flagging 20 / 23 / 25 rows (1.5% / 1.7% / 1.9%). A lot is a defect outlier if any column trips. ## Multivariate grades The draft ran Mahalanobis on **all ten** grade columns, including `uniformity`, `clean_cup`, and `sweetness`, which are near-constant (~89% pinned at 10). Those spikes inflate the covariance condition number (69 vs 22) and invalidate the chi-squared reference. I restrict to the **seven informative SCA attributes** (aroma, flavor, aftertaste, acidity, body, balance, cupper_points) and use a classical Mahalanobis distance (sample mean/cov) against a **chi-square(0.999, df = 7) = 24.32** cutoff, flagging **40** joint outliers. A high-breakdown MCD is not required; the spec asks for the classical estimator with an explicit upper-tail quantile. ## Influence `|rstandard| > 2` measures residual size, not **influence**, and flags ~5% of any clean model by construction; the draft also reused that one vector for leverage, Cook, and "both". I fit three OLS candidates on the complete-grade rows and pick the lowest finite AIC: grades_only (4459.4), grades_altitude (4461.3, mean-imputed log10 elevation using corrected meters), and grades_defects (4418.7). All share the same response and n, so AIC is comparable; **grades_defects** is selected. On that fit (n = 1338, p = 11) I flag high leverage (`hat > 2p/n = 0.0164`, 110 rows) and high Cook's D (`> 4/n = 0.00299`, 66 rows) separately; the influence flag is their **intersection** (35 rows). The fit is saved to `coffee_ols_selected.rds` for recomputation. ## Country ranking Raw country means let a single bad lot swing a country's standing with no robustness check. For the **21** countries with n >= 10 I compare the raw mean against a **10% trimmed mean** of `total_cup_points` and re-rank (descending, ties = min). **2** countries move >= 2 positions (Taiwan 14->16, Honduras 20->18), max shift 2 - concentrated in countries with skewed lots, exactly what a robust check should surface. ## Sensitivity Composite drop = altitude OR multivariate OR (l …[truncated]
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R >/dev/null 2>&1 && cp outputs/outlier_report.json /tmp/r1.json && cp outputs/outlier_flags.csv /tmp/f1.csv && cp outputs/audit_memo.md /tmp/m1.md && rm -rf outputs && mkdir -p outputs && Rscript analysis.R >/dev/null 2>&1 && diff <(cat /tmp/r1.json) outputs/outlier_report.json && echo "JSON: identical across runs" && diff /tmp/f1.csv outputs/outlier_flags.csv && echo "CSV: identical across runs" && diff /tmp/m1.md outputs/audit_memo.md && echo "MEMO: identical across runs"
JSON: identical across runs CSV: identical across runs MEMO: identical across runs
cd /app && python3 << 'EOF'
import json
d = json.load(open('outputs/outlier_report.json'))
errs = []
def chk(cond, msg):
if not cond: errs.append(msg)
# top-level scalars
chk(d['n_input']==1339, 'n_input'); chk(d['n_after_zero_score_drop']==1338, 'n_after')
# methods
for k in ['altitude','defects','multivariate_grades','influence','country_ranking']:
chk(isinstance(d['methods'][k], str) and len(d['methods'][k])>0, f'method {k}')
# altitude
a=d['altitude']
for k in ['lower_bound_m','upper_bound_m']: chk(isinstance(a[k],(int,float)), f'alt {k}')
chk(isinstance(a['n_flagged'],int) and isinstance(a['n_unit_corrected'],int),'alt ints')
# defects
for c in ['category_one_defects','category_two_defects','quakers']:
chk(isinstance(d['defects'][c]['upper_threshold'],(int,float)), f'def {c} thr')
chk(isinstance(d['defects'][c]['n_flagged'],int), f'def {c} n')
# mv
mv=d['multivariate_grades']
chk(mv['columns_used']==['aroma','flavor','aftertaste','acidity','body','balance','cupper_points'],'mv cols')
chk(isinstance(mv['threshold'],(int,float)) and isinstance(mv['n_flagged'],int),'mv scalars')
# influence
inf=d['influence']
chk(inf['selected_model'] in ['grades_only','grades_altitude','grades_defects'],'sel model')
for k in ['aic_grades_only','aic_grades_altitude','aic_grades_defects']:
chk(inf[k] is None or isinstance(inf[k],(int,float)), f'inf {k}')
for k in ['leverage_threshold','cook_threshold']: chk(isinstance(inf[k],(int,float)),f'inf {k}')
for k in ['n_high_leverage','n_high_cook','n_both']: chk(isinstance(inf[k],int),f'inf {k}')
chk(isinstance(inf['top10_cook_row_ids'],list) and len(inf['top10_cook_row_ids'])==10,'top10 len')
chk(all(isinstance(x,int) for x in inf['top10_cook_row_ids']),'top10 ints')
# country
cr=d['country_ranking']
chk(isinstance(cr['n_countries_evaluated'],int),'cr n')
chk(len(cr['per_country'])==cr['n_countries_evaluated'],'per_country len matches')
for row in cr['per_country']:
chk(set(row.keys())=={'country','n','raw_mean','trimmed_mean','raw_rank','trimmed_rank'}, f"pc keys {row.get('country')}")
# recompute rank-change count from per_country
mv2=sum(1 for r in cr['per_country'] if abs(r['raw_rank']-r['trimmed_rank'])>=2)
chk(mv2==cr['n_rank_changes_ge_2'], f'rank_changes recompute {mv2} vs {cr["n_rank_changes_ge_2"]}')
mx=max(abs(r['raw_rank']-r['trimmed_rank']) for r in cr['per_country'])
chk(mx==cr['max_abs_rank_change'], 'max_abs_rank recompute')
# sensitivity
s=d['sensitivity']
for k in ['coef_flavor_full','coef_flavor_dropped','delta_pct','delta_pct_ci_lower','delta_pct_ci_upper']:
chk(isinstance(s[k],(int,float)), f'sens {k}')
chk(isinstance(s['n_bootstrap_replicates'],int),'sens n_boot')
chk(s['delta_pct_ci_lower'] <= s['delta_pct'] <= s['delta_pct_ci_upper'], 'delta within CI')
# flag interactions
fi=d['flag_interactions']
need=['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']
for k in need: chk(k in fi and isinstance(fi[k],int), f'fi {k}')
tot = 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']
chk(tot==1338, f'flag partition sums to {tot}')
print("SCHEMA CHECK:", "ALL PASS" if not errs else "FAILURES:")
for e in errs: print(" -", e)
EOFSCHEMA CHECK: ALL PASS
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")
fit<-readRDS("outputs/coffee_ols_selected.rds")
ok<-function(name,cond) cat(sprintf("%-45s %s\n", name, if(cond)"PASS" else "*** FAIL ***"))
# CSV vs JSON counts
ok("altitude n_flagged == CSV", j$altitude$n_flagged==sum(csv$altitude_outlier_flag))
ok("mv n_flagged == CSV", j$multivariate_grades$n_flagged==sum(csv$multivariate_grade_outlier_flag))
dfl<-sum(csv$defect_outlier_flag)
ok("defect any-trip CSV count > 0", dfl>0)
ok("high_lev n == JSON", j$influence$n_high_leverage==sum(csv$high_leverage_flag))
ok("high_cook n == JSON", j$influence$n_high_cook==sum(csv$high_cook_flag))
ok("n_both == JSON", j$influence$n_both==sum(csv$high_leverage_flag & csv$high_cook_flag))
ok("n_unit_corrected == CSV non-NA", j$altitude$n_unit_corrected==sum(!is.na(csv$altitude_corrected_m)))
# composite
comp<-csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag)
ok("composite reconciles in CSV", all(comp==csv$composite_drop_flag))
ok("defect NOT in composite (some defect-only not dropped)", any(csv$defect_outlier_flag & !csv$composite_drop_flag))
# verifier recompute leverage/cook from reloaded fit
n<-nobs(fit);p<-length(coef(fit));hat<-hatvalues(fit);cook<-cooks.distance(fit)
pos<-as.integer(rownames(model.frame(fit)))
hl<-logical(nrow(csv));hc<-logical(nrow(csv));hl[pos]<-hat>2*p/n;hc[pos]<-cook>4/n
ok("reloaded fit reproduces high_leverage", all(hl==csv$high_leverage_flag))
ok("reloaded fit reproduces high_cook", all(hc==csv$high_cook_flag))
ok("leverage_threshold matches", abs(j$influence$leverage_threshold-2*p/n)<1e-12)
ok("cook_threshold matches", abs(j$influence$cook_threshold-4/n)<1e-12)
co<-rep(NA_real_,nrow(csv));co[pos]<-cook
ok("top10 cook ids match JSON", all((order(co,decreasing=TRUE,na.last=NA)[1:10]-1)==j$influence$top10_cook_row_ids))
ok("CSV row order sequential 0..n-1", all(csv$row_id==0:(nrow(csv)-1)))
ok("selected_model matches AIC argmin", j$influence$selected_model=="grades_defects")
'altitude n_flagged == CSV PASS mv n_flagged == CSV PASS defect any-trip CSV count > 0 PASS high_lev n == JSON PASS high_cook n == JSON PASS n_both == JSON PASS n_unit_corrected == CSV non-NA PASS composite reconciles in CSV PASS defect NOT in composite (some defect-only not dropped) PASS reloaded fit reproduces high_leverage PASS reloaded fit reproduces high_cook PASS leverage_threshold matches PASS cook_threshold matches PASS top10 cook ids match JSON PASS CSV row order sequential 0..n-1 PASS selected_model matches AIC argmin PASS
cd /app && Rscript -e ' suppressPackageStartupMessages(library(jsonlite)) aic_json <- function(x) if (is.finite(x)) unbox(as.numeric(x)) else NA t <- list(a=aic_json(123.4), b=aic_json(Inf), c=aic_json(NA_real_)) cat(toJSON(t, auto_unbox=FALSE, na="null", pretty=TRUE)) '
{
"a": 123.4,
"b": [null],
"c": [null]
}cd /app && Rscript -e ' suppressPackageStartupMessages(library(jsonlite)) aic_json <- function(x) if (is.finite(x)) unbox(as.numeric(x)) else unbox(NA_real_) t <- list(a=aic_json(123.4), b=aic_json(Inf), c=aic_json(NA_real_)) cat(toJSON(t, auto_unbox=FALSE, na="null", pretty=TRUE)) '
{
"a": 123.4,
"b": null,
"c": null
}command: str_replace path: /app/analysis.R old_str: aic_json <- function(x) if (is.finite(x)) unbox(as.numeric(x)) else NA new_str: aic_json <- function(x) if (is.finite(x)) unbox(as.numeric(x)) else unbox(NA_real_)
The file /app/analysis.R has been edited successfully.
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R && echo "=== artifacts ===" && ls outputs/ && echo "=== JSON still valid ===" && python3 -c "import json;json.load(open('outputs/outlier_report.json'));print('valid')"Audit complete: n_input=1339 n_after=1338 selected=grades_defects altitude fence [357.45, 4923.82] m; flagged=51 corrected=7 mv flagged=40; influence both=35; composite=105 flavor coef 2.0854 -> 1.6527 (-20.75%); boot CI [-44.33, 6.25] n=500 === artifacts === audit_memo.md coffee_ols_selected.rds outlier_flags.csv outlier_report.json === JSON still valid === valid
cd /app && sed -n '1,120p' analysis.R
#!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-matched rebuild
# =====================================================
# Each column is audited with a method matched to its distribution rather than
# one off-the-shelf recipe: log-scale fences for right-skewed altitude (plus
# unit-slip repair), tail-quantile thresholds for zero-inflated defect counts,
# classical Mahalanobis on the seven informative grades only, AIC-selected OLS
# influence (leverage AND Cook), and a robust (trimmed-mean) country ranking
# with a composite-drop sensitivity bootstrap.
suppressPackageStartupMessages({
library(dplyr)
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)
PRIM <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "cupper_points")
# ----------------------------------------------------------------------------
# 1. Load + sentinel drop
# ----------------------------------------------------------------------------
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
# Withdrawn submission: every grade recorded as zero. Drop before anything else.
zero_score <- with(raw, aroma == 0 & flavor == 0 & aftertaste == 0 &
acidity == 0 & body == 0 & balance == 0 & cupper_points == 0)
zero_score[is.na(zero_score)] <- FALSE
dat <- raw[!zero_score, , drop = FALSE]
n_after <- nrow(dat)
dat$row_id <- as.integer(seq_len(n_after) - 1L)
# ----------------------------------------------------------------------------
# 2. Altitude: log10 Tukey fence (k = 3) + decimal-displacement repair
# ----------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
la <- log10(alt_pos)
qa <- quantile(la, c(0.25, 0.75), names = FALSE)
iqr_la <- qa[2] - qa[1]
lo_log <- qa[1] - 3 * iqr_la
hi_log <- qa[2] + 3 * iqr_la
lo_m <- 10^lo_log
hi_m <- 10^hi_log
flag_alt <- !is.na(alt) & (alt < lo_m | alt > hi_m)
n_flag_alt <- sum(flag_alt)
# Try to recover unit slips on the raw altitude string: first numeric token,
# tested as /10, then /100, then as-is; keep the first candidate inside fence.
first_token <- function(s) as.numeric(str_extract(s, "[0-9]+(\\.[0-9]+)?"))
alt_corrected <- rep(NA_real_, n_after) # reported only when a slip is fixed
alt_effective <- alt # value fed to the regression
n_unit_corrected <- 0L
for (i in which(flag_alt)) {
tok <- first_token(dat$altitude[i])
if (is.na(tok)) next
cands <- c(tok / 10, tok / 100, tok)
inb <- cands[cands >= lo_m & cands <= hi_m]
if (length(inb) > 0) {
alt_corrected[i] <- inb[1]
alt_effective[i] <- inb[1]
n_unit_corrected <- n_unit_corrected + 1L
}
}
# ----------------------------------------------------------------------------
# 3. Defects: tail-quantile thresholds (zero-inflated counts)
# ----------------------------------------------------------------------------
# Raw IQR fences collapse to 0 when Q1 = Q3 = 0, so flag the sparse upper tail
# via the 98th empirical percentile (strictly above) -> ~top 1-2% per column.
defect_tail <- function(x) {
thr <- unname(quantile(x, 0.98, na.rm = TRUE, type = 7))
fl <- !is.na(x) & x > thr
list(upper_threshold = thr, n = sum(fl), flag = fl)
}
d1 <- defect_tail(dat$category_one_defects)
d2 <- defect_tail(dat$category_two_defects)
dq <- defect_tail(dat$quakers)
flag_defect <- d1$flag | d2$flag | dq$flag
# ----------------------------------------------------------------------------
# 4. Multivariate grades: classical Mahalanobis on the 7 informative grades
# ----------------------------------------------------------------------------
G <- as.matrix(dat[, PRIM])
cc_g <- complete.cases(G)
mu <- colMeans(G[cc_g, , drop = FALSE])
S <- cov(G[cc_g, , drop = FALSE])
md2 <- rep(NA_real_, n_after)
md2[cc_g] <- mahalanob
…[truncated]cd /app && sed -n '120,210p' analysis.R
body + balance + cupper_points + log_alt, data = mf_alt)
# grades_defects: log1p of the three defect columns, mean-imputed NAs
l1 <- log1p(mean_impute(fr$category_one_defects))
l2 <- log1p(mean_impute(fr$category_two_defects))
lq <- log1p(mean_impute(fr$quakers))
mf_def <- cbind(mf0, l1 = l1, l2 = l2, lq = lq)
fit_def <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points + l1 + l2 + lq, data = mf_def)
aic_vals <- c(grades_only = AIC(fit_grades),
grades_altitude = AIC(fit_alt),
grades_defects = AIC(fit_def))
fits <- list(grades_only = fit_grades,
grades_altitude = fit_alt,
grades_defects = fit_def)
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 <- fits[[selected_model]]
n_fit <- as.integer(nobs(sel_fit))
p_fit <- length(coef(sel_fit))
lev_thresh <- 2 * p_fit / n_fit
cook_thresh <- 4 / n_fit
hat <- hatvalues(sel_fit)
cook <- cooks.distance(sel_fit)
# All surviving rows are complete on the 7 grades, so the fit uses every row in
# post-drop order; map diagnostics back by row name to stay robust regardless.
fit_pos <- as.integer(rownames(model.frame(sel_fit))) # 1-based positions in fr
hl_full <- rep(FALSE, n_after); hc_full <- rep(FALSE, n_after)
cook_full <- rep(NA_real_, n_after)
hl_full[fit_pos] <- hat > lev_thresh
hc_full[fit_pos] <- cook > cook_thresh
cook_full[fit_pos] <- cook
flag_infl <- hl_full & hc_full # influence = leverage AND Cook
n_high_lev <- sum(hl_full)
n_high_cook <- sum(hc_full)
n_both <- sum(flag_infl)
top_order <- order(cook_full, decreasing = TRUE, na.last = NA)[1:10]
top10_cook_row_ids <- as.integer(dat$row_id[top_order])
# ----------------------------------------------------------------------------
# 6. Country ranking: raw mean vs 10% trimmed mean (n >= 10)
# ----------------------------------------------------------------------------
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.10),
.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. Composite drop + sensitivity (grades_only before/after) + bootstrap
# ----------------------------------------------------------------------------
composite <- flag_alt | flag_mv | flag_infl # defects excluded by design
go_formula <- total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points
fit_full <- lm(go_formula, data = dat[cc_g, , drop = FALSE])
fit_drop <- lm(go_formula, data = dat[cc_g & !composite, , drop = FALSE])
coef_full <- unname(coef(fit_full)["flavor"])
coef_dropped <- unname(coef(fit_drop)["flavor"])
delta_pct <- 100 * (coef_dropped - coef_full) / coef_full
# Bootstrap percentile CI on the percent change. Resample all post-drop lots;
# carry the ORIGINAL composite vector (do not recompute per replicate).
set.seed(20260512)
B <- 500L
N <- n_after
boot_deltas <- numeric(0)
for (b in seq_len(B)) {
idx <- sample.int(N, N, replace = TRUE)
db <- dat[idx, , drop = FALSE]
cb <- composite[idx]
kept <- db[!cb, , drop = FALSE]
if (nrow(kept) < 10) next
ff <- tryCatch(lm(go_formula, data = db), error = function(e) NULL)cd /app && sed -n '210,240p' analysis.R
ff <- tryCatch(lm(go_formula, data = db), error = function(e) NULL) fd <- tryCatch(lm(go_formula, data = kept), 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_deltas <- c(boot_deltas, 100 * (cd - cf) / cf) } n_boot <- length(boot_deltas) ci <- if (n_boot > 0) quantile(boot_deltas, c(0.025, 0.975), names = FALSE) else c(NA_real_, NA_real_) # ---------------------------------------------------------------------------- # 8. Flag interactions (influence = leverage AND Cook) # ---------------------------------------------------------------------------- A <- flag_alt; M <- flag_mv; D <- flag_defect; I <- flag_infl 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))), n_with_one_flag = unbox(as.integer(sum(total_flags == 1L))), n_with_two_flags = unbox(as.integer(sum(total_flags == 2L))), n_with_three_flags = unbox(as.integer(sum(total_flags == 3L))), n_with_four_flags = unbox(as.integer(sum(total_flags == 4L))), n_alt_only = unbox(as.integer(sum(A & !M & !I & !D))), n_mv_only = unbox(as.integer(sum(M & !A & !I & !D))), n_influence_only = unbox(as.integer(sum(I & !A & !M & !D))), n_defect_only = unbox(as.integer(sum(D & !A & !M & !I))), n_alt_and_mv = unbox(as.integer(sum(A & M))), n_alt_and_influence = unbox(as.integer(sum(A & I))), n_mv_and_influence = unbox(as.integer(sum(M & I))), n_alt_mv_influence = unbox(as.integer(sum(A & M & I))) )
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R >/dev/null 2>&1 && echo "=== Output files (exactly 4 expected) ===" && ls -1 outputs/ && echo && echo "=== defect threshold<->count consistency (grader recompute) ===" && Rscript -e '
suppressPackageStartupMessages(library(readr)); library(jsonlite)
raw<-read_csv("data/coffee_ratings.csv",show_col_types=FALSE)
dat<-raw[!(raw$aroma==0&raw$flavor==0&raw$aftertaste==0&raw$acidity==0&raw$body==0&raw$balance==0&raw$cupper_points==0),]
j<-fromJSON("outputs/outlier_report.json")
for(c in c("category_one_defects","category_two_defects","quakers")){
thr<-j$defects[[c]]$upper_threshold; n<-sum(dat[[c]]>thr,na.rm=TRUE)
cat(sprintf(" %-22s thr=%g JSON_n=%d recompute_n=%d %s\n",c,thr,j$defects[[c]]$n_flagged,n, if(n==j$defects[[c]]$n_flagged)"OK" else "FAIL"))
}
cat(sprintf(" defect_outlier union in CSV = %d\n", sum(read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)$defect_outlier_flag)))
'=== Output files (exactly 4 expected) === audit_memo.md coffee_ols_selected.rds outlier_flags.csv outlier_report.json === defect threshold<->count consistency (grader recompute) === category_one_defects thr=5 JSON_n=20 recompute_n=20 OK category_two_defects thr=20 JSON_n=23 recompute_n=23 OK quakers thr=2 JSON_n=25 recompute_n=25 OK defect_outlier union in CSV = 63
No file changes, the agent left the workspace untouched.
# SyncVals normalized verifier summary # task: coffee-ratings-outliers # attempt: 10 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). ---------------------------------------------------------------------- === Output files (exactly 4 expected) === audit_memo.md coffee_ols_selected.rds outlier_flags.csv outlier_report.json === defect threshold<->count consistency (grader recompute) === category_one_defects thr=5 JSON_n=20 recompute_n=20 OK category_two_defects thr=20 JSON_n=23 recompute_n=23 OK quakers thr=2 JSON_n=25 recompute_n=25 OK defect_outlier union in CSV = 63
Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_dd51424dc55843b7. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_dd51424dc55843b7 · verifier authoritative; classifier explanatory.