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 2>/dev/null
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 "---COLNAMES---" && 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
---COLNAMES---
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 "---PKGS---" && Rscript -e 'for (p in c("dplyr","tidyr","readr","stringr","jsonlite","purrr")) cat(p, as.character(packageVersion(p)), "\n")' 2>&1R version 4.4.2 (2024-10-31) -- "Pile of Leaves" ---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
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")
# the all-zero row
allz <- raw %>% filter(if_all(all_of(grades), ~ .x == 0))
cat("rows all grades zero:", nrow(allz), "\n")
cat("total_cup_points==0 rows:", sum(raw$total_cup_points==0, na.rm=TRUE), "\n")
print(raw %>% filter(total_cup_points==0) %>% select(total_cup_points, all_of(grades)))
cat("\n--- summary of grades ---\n")
print(summary(raw[grades]))
' 2>&1n_input: 1339
rows all grades zero: 1
total_cup_points==0 rows: 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>
--- summary of grades ---
aroma flavor aftertaste acidity body
Min. :0.000 Min. :0.00 Min. :0.000 Min. :0.000 Min. :0.000
1st Qu.:7.420 1st Qu.:7.33 1st Qu.:7.250 1st Qu.:7.330 1st Qu.:7.330
Median :7.580 Median :7.58 Median :7.420 Median :7.580 Median :7.500
Mean :7.567 Mean :7.52 Mean :7.401 Mean :7.536 Mean :7.517
3rd Qu.:7.750 3rd Qu.:7.75 3rd Qu.:7.580 3rd Qu.:7.750 3rd Qu.:7.670
Max. :8.750 Max. :8.83 Max. :8.670 Max. :8.750 Max. :8.580
balance uniformity clean_cup sweetness
Min. :0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.:7.330 1st Qu.:10.000 1st Qu.:10.000 1st Qu.:10.000
Median :7.500 Median :10.000 Median :10.000 Median :10.000
Mean :7.518 Mean : 9.835 Mean : 9.835 Mean : 9.857
3rd Qu.:7.750 3rd Qu.:10.000 3rd Qu.:10.000 3rd Qu.:10.000
Max. :8.750 Max. :10.000 Max. :10.000 Max. :10.000
cupper_points
Min. : 0.000
1st Qu.: 7.250
Median : 7.500
Mean : 7.503
3rd Qu.: 7.750
Max. :10.000cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
for (c in c("category_one_defects","category_two_defects","quakers")) {
x <- raw[[c]]
cat("===", c, "=== NA:", sum(is.na(x)), " zeros:", sum(x==0,na.rm=TRUE), "/", length(x), "\n")
print(quantile(x, c(.5,.9,.95,.97,.98,.99,1), na.rm=TRUE))
cat("IQR Q1/Q3:", quantile(x,.25,na.rm=TRUE), quantile(x,.75,na.rm=TRUE), "\n\n")
}
' 2>&1=== category_one_defects === NA: 0 zeros: 1137 / 1339 50% 90% 95% 97% 98% 99% 100% 0.00 1.00 2.00 4.00 5.00 7.62 63.00 IQR Q1/Q3: 0 0 === category_two_defects === NA: 0 zeros: 373 / 1339 50% 90% 95% 97% 98% 99% 100% 2.00 9.00 13.00 17.00 20.00 26.62 55.00 IQR Q1/Q3: 0 4 === quakers === NA: 1 zeros: 1244 / 1339 50% 90% 95% 97% 98% 99% 100% 0 0 1 2 2 5 11 IQR 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)
for (c in c("category_one_defects","category_two_defects","quakers")) {
x <- dat[[c]]; x <- x[!is.na(x)]
nz <- x[x>0]
cat("===", c, "=== n:", length(x), " frac>0:", round(mean(x>0),3), "\n")
# Poisson upper tail
pt <- qpois(0.99, mean(x))
# quantile thresholds
q98 <- quantile(x, .98); q99 <- quantile(x, .99)
# Tukey on nonzero
qn <- quantile(nz, c(.25,.75)); fen <- qn[[2]] + 1.5*(qn[[2]]-qn[[1]]); fen3 <- qn[[2]] + 3*(qn[[2]]-qn[[1]])
cat(sprintf(" qpois99=%g (flag>%g -> n=%d, %.1f%%)\n", pt, pt, sum(x>pt), 100*mean(x>pt)))
cat(sprintf(" q98=%g (n=%d, %.1f%%) q99=%g (n=%d, %.1f%%)\n", q98, sum(x>q98),100*mean(x>q98), q99, sum(x>q99),100*mean(x>q99)))
cat(sprintf(" Tukey-nonzero 1.5=%g (n=%d, %.1f%%) 3.0=%g (n=%d, %.1f%%)\n", fen, sum(x>fen),100*mean(x>fen), fen3, sum(x>fen3),100*mean(x>fen3)))
}
' 2>&1=== category_one_defects === n: 1338 frac>0: 0.151 qpois99=3 (flag>3 -> n=45, 3.4%) q98=5 (n=20, 1.5%) q99=7.63 (n=14, 1.0%) Tukey-nonzero 1.5=6 (n=17, 1.3%) 3.0=9 (n=12, 0.9%) === category_two_defects === n: 1338 frac>0: 0.721 qpois99=9 (flag>9 -> n=115, 8.6%) q98=20 (n=23, 1.7%) q99=26.63 (n=14, 1.0%) Tukey-nonzero 1.5=12 (n=73, 5.5%) 3.0=18 (n=34, 2.5%) === quakers === n: 1337 frac>0: 0.07 qpois99=2 (flag>2 -> n=25, 1.9%) q98=2 (n=25, 1.9%) q99=5 (n=10, 0.7%) Tukey-nonzero 1.5=6 (n=6, 0.4%) 3.0=9 (n=1, 0.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)
alt <- dat$altitude_mean_meters
cat("NA:", sum(is.na(alt)), " <=0:", sum(alt<=0, na.rm=TRUE), " >0:", sum(alt>0,na.rm=TRUE), "\n")
ap <- alt[!is.na(alt) & alt>0]
cat("range meters:", range(ap), "\n")
print(quantile(ap, c(0,.25,.5,.75,.9,.99,1)))
cat("\n--- log10 fence k=3 ---\n")
la <- log10(ap)
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,IQR:", q[[1]], q[[2]], iqr, "\n")
cat("log fence:", lo, hi, " -> meters:", 10^lo, 10^hi, "\n")
flag <- !is.na(alt) & alt>0 & (log10(alt) < lo | log10(alt) > hi)
# also rows with NA or <=0? flag only out of fence among positive
cat("n flagged (positive, outside):", sum(flag), "\n")
cat("rows alt<=0:", sum(!is.na(alt) & alt<=0), "\n")
cat("\nHigh altitudes flagged:\n")
print(dat %>% filter(flag) %>% select(altitude, altitude_mean_meters) %>% arrange(desc(altitude_mean_meters)) %>% head(40), n=40)
' 2>&1NA: 230 <=0: 0 >0: 1108
range meters: 1 190164
0% 25% 50% 75% 90% 99% 100%
1.00 1100.00 1310.64 1600.00 1799.30 3280.00 190164.00
--- log10 fence k=3 ---
log10 Q1,Q3,IQR: 3.041393 3.20412 0.1627273
log fence: 2.553211 3.692302 -> meters: 357.4463 4923.817
n flagged (positive, outside): 51
rows alt<=0: 0
High altitudes flagged:
# A tibble: 40 × 2
altitude altitude_mean_meters
<chr> <dbl>
1 190164 190164
2 1901.64 190164
3 1100.00 mosl 110000
4 11000 metros 11000
5 350 350
6 350 350
7 ~350m 350
8 300 m 300
9 280 280
10 approx. 1'500 metres a.s.l. 250.
11 approx.1'500 metres a.s.l. 250.
12 250 m 250
13 250 m 250
14 200 m 200
15 200 m 200
16 200 200
17 200 m 200
18 200 200
19 1.85 185
20 160-200m 180
21 175 175
22 170 170
23 16.8 168
24 1.65 165
25 518 158.
26 518 158.
27 518 158.
28 150 m 150
29 150 m 150
30 1.25 125
31 110公尺 110
32 100 100
33 50 50
34 40 40
35 1.3 13
36 1.3 13
37 1.2 12
38 1.2 12
39 12oo 12
40 -1 1cd /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]; la <- log10(ap)
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(alt) & alt>0 & (alt < lo | alt > hi)
first_num <- function(s){ m <- str_extract(s, "[0-9]+(\\.[0-9]+)?"); as.numeric(m) }
corr <- function(raw_str){
tok <- first_num(raw_str)
if (is.na(tok)) return(NA_real_)
for (cand in c(tok/10, tok/100, tok)) if (!is.na(cand) && cand>=lo && cand<=hi) return(cand)
NA_real_
}
fl <- which(flag)
res <- sapply(dat$altitude[fl], corr)
cat(sprintf("fence lo=%.3f hi=%.3f n_flag=%d n_corrected=%d\n", lo, hi, length(fl), sum(!is.na(res))))
df <- tibble(altitude=dat$altitude[fl], mean_m=alt[fl], tok=sapply(dat$altitude[fl], first_num), corrected=res)
print(df, n=60)
' 2>&1fence lo=357.446 hi=4923.817 n_flag=51 n_corrected=7 # A tibble: 51 × 4 altitude mean_m tok 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 <- mahalanobis(G[cc,], mu, S)
th <- qchisq(0.999, df=7)
cat("threshold:", th, " n_flagged:", sum(md2>th), "\n")
# model AICs
dat$row_id <- seq_len(nrow(dat))-1L
fit_df <- dat[cc,]
f1 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points, data=fit_df)
cat("AIC grades_only:", AIC(f1), " n:", nobs(f1), " p:", length(coef(f1)), "\n")
cat("R2:", summary(f1)$r.squared, "\n")
# residual sd
cat("sigma:", summary(f1)$sigma, "\n")
' 2>&1complete cases on 7 grades: 1338 of 1338 threshold: 24.32189 n_flagged: 40 AIC grades_only: 4459.362 n: 1338 p: 8 R2: 0.7756335 sigma: 1.27604
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 corrected
alt <- dat$altitude_mean_meters
ap <- alt[!is.na(alt)&alt>0]; la <- log10(ap); 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(alt)&alt>0&(alt<lo|alt>hi)
first_num<-function(s){as.numeric(str_extract(s,"[0-9]+(\\.[0-9]+)?"))}
corr<-rep(NA_real_,nrow(dat))
for(i in which(flag)){tok<-first_num(dat$altitude[i]); if(!is.na(tok)) for(cand in c(tok/10,tok/100,tok)) if(cand>=lo&cand<=hi){corr[i]<-cand;break}}
alt_reg <- ifelse(!is.na(corr), corr, alt)
G <- as.matrix(dat[,g7]); cc <- complete.cases(G)
fit_df <- dat[cc,]
# model 1
f1 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points, data=fit_df)
# model 2: + log10 elevation, mean-impute nonfinite
logelev <- log10(alt_reg[cc]); logelev[!is.finite(logelev)] <- NA
logelev[is.na(logelev)] <- mean(logelev, na.rm=TRUE)
fit_df2 <- fit_df; fit_df2$log_elev <- logelev
f2 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+log_elev, data=fit_df2)
# model 3: + log1p defects, mean-impute NA
d1<-dat$category_one_defects[cc]; d2<-dat$category_two_defects[cc]; dq<-dat$quakers[cc]
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
fit_df3<-fit_df; fit_df3$ld1<-log1p(imp(d1)); fit_df3$ld2<-log1p(imp(d2)); fit_df3$ldq<-log1p(imp(dq))
f3 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+ld1+ld2+ldq, data=fit_df3)
cat(sprintf("AIC grades_only=%.4f\nAIC grades_altitude=%.4f\nAIC grades_defects=%.4f\n", AIC(f1),AIC(f2),AIC(f3)))
cat("nobs:", nobs(f1), nobs(f2), nobs(f3), "\n")
' 2>&1AIC grades_only=4459.3615 AIC grades_altitude=4461.2923 AIC grades_defects=4418.6801 nobs: 1338 1338 1338
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); fit_df <- dat[cc,]
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
fit_df$ld1<-log1p(imp(dat$category_one_defects[cc])); fit_df$ld2<-log1p(imp(dat$category_two_defects[cc])); fit_df$ldq<-log1p(imp(dat$quakers[cc]))
f3 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+ld1+ld2+ldq, data=fit_df)
n<-nobs(f3); p<-length(coef(f3))
lev<-hatvalues(f3); ck<-cooks.distance(f3)
lev_th<-2*p/n; ck_th<-4/n
cat("n:",n," p:",p," lev_th:",lev_th," cook_th:",ck_th,"\n")
cat("n_high_lev:",sum(lev>lev_th)," n_high_cook:",sum(ck>ck_th)," n_both:",sum(lev>lev_th & ck>ck_th),"\n")
# top10 cook row ids (0-indexed, positions among post-drop; cc all TRUE so positions=1:n)
ord<-order(-ck)[1:10]
cat("top10 cook 0-idx row_ids:", paste((which(cc)[ord]-1L), collapse=","), "\n")
cat("top10 cook values:", paste(round(ck[ord],4),collapse=","), "\n")
' 2>&1n: 1338 p: 11 lev_th: 0.01644245 cook_th: 0.002989537 n_high_lev: 110 n_high_cook: 66 n_both: 35 top10 cook 0-idx row_ids: 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"),
d=abs(raw_rank-trimmed_rank))
cat("=== trim",tr,"=== n_countries:",nrow(ct)," n_changes>=2:",sum(ct$d>=2)," max:",max(ct$d),"\n")
}
# show detail for trim 0.1
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"), d=abs(raw_rank-trimmed_rank)) %>% arrange(raw_rank)
print(ct, n=40)
' 2>&1=== trim 0.1 === n_countries: 21 n_changes>=2: 2 max: 2 === trim 0.2 === n_countries: 21 n_changes>=2: 2 max: 2 # A tibble: 21 × 7 country n raw_mean trimmed_mean raw_rank trimmed_rank d <chr> <int> <dbl> <dbl> <int> <int> <int> 1 Ethiopia 44 85.5 85.5 1 1 0 2 United States 10 84.4 84.9 2 2 0 3 Kenya 25 84.3 84.5 3 3 0 4 Uganda 36 83.5 83.4 4 4 0 5 Colombia 183 83.1 83.2 5 5 0 6 El Salvador 21 83.1 83.1 6 6 0 7 China 16 82.9 83.0 7 8 1 8 Costa Rica 51 82.8 83.0 8 7 1 9 Thailand 32 82.6 82.6 9 10 1 10 Indonesia 20 82.6 82.8 10 9 1 11 Peru 10 82.5 82.4 11 12 1 12 Brazil 132 82.4 82.5 12 11 1 13 Tanzania, United Rep… 40 82.4 82.3 13 13 0 14 Taiwan 75 82.0 82.0 14 16 2 15 Guatemala 181 81.8 82.2 15 14 1 16 United States (Hawai… 73 81.8 82.1 16 15 1 17 Malawi 11 81.7 81.7 17 17 0 18 India 14 81.1 81.4 18 19 1 19 Mexico 236 80.9 81.2 19 20 1 20 Honduras 52 80.9 81.5 20 18 2 21 Nicaragua 26 80.5 80.9 21 21 0
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")
G <- as.matrix(dat[,g7]); cc <- complete.cases(G)
# altitude flag
alt<-dat$altitude_mean_meters; ap<-alt[!is.na(alt)&alt>0]; la<-log10(ap); 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(alt)&alt>0&(alt<lo|alt>hi)
# mv flag
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(0.999,7)
# selected model infl
fit_df<-dat[cc,]; imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
fit_df$ld1<-log1p(imp(dat$category_one_defects[cc])); fit_df$ld2<-log1p(imp(dat$category_two_defects[cc])); fit_df$ldq<-log1p(imp(dat$quakers[cc]))
f3<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+ld1+ld2+ldq,data=fit_df)
n3<-nobs(f3);p3<-length(coef(f3)); lev<-rep(FALSE,nrow(dat));ckf<-rep(FALSE,nrow(dat))
lev[cc]<-hatvalues(f3)>2*p3/n3; ckf[cc]<-cooks.distance(f3)>4/n3
composite <- flag_alt | flag_mv | (lev & ckf)
cat("composite n:", sum(composite), " (alt",sum(flag_alt),"mv",sum(flag_mv),"both",sum(lev&ckf),")\n")
# grades_only full vs dropped
f_full<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat[cc,])
keep <- cc & !composite
f_drop<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat[keep,])
cf<-unname(coef(f_full)["flavor"]); cd<-unname(coef(f_drop)["flavor"])
dp<-100*(cd-cf)/cf
cat(sprintf("coef_full=%.6f coef_drop=%.6f delta_pct=%.6f\n",cf,cd,dp))
# bootstrap
set.seed(20260512); B<-500
df_all <- dat[cc,]; comp_all <- composite[cc]; nboot<-nrow(df_all)
deltas<-numeric(0)
for(b in 1:B){
idx<-sample.int(nboot,nboot,replace=TRUE)
sub<-df_all[idx,]; comp_b<-comp_all[idx]
dsub<-sub[!comp_b,]
if(nrow(dsub)<10) next
ff<-tryCatch(lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=sub),error=function(e)NULL)
fd<-tryCatch(lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dsub),error=function(e)NULL)
if(is.null(ff)||is.null(fd)) next
cfb<-unname(coef(ff)["flavor"]); cdb<-unname(coef(fd)["flavor"])
if(is.na(cfb)||is.na(cdb)||cfb==0) next
deltas<-c(deltas,100*(cdb-cfb)/cfb)
}
ci<-quantile(deltas,c(.025,.975))
cat(sprintf("n_boot_success=%d CI=[%.4f, %.4f]\n",length(deltas),ci[1],ci[2]))
' 2>&1composite n: 105 (alt 51 mv 40 both 35 ) coef_full=2.085398 coef_drop=1.652728 delta_pct=-20.747640 n_boot_success=500 CI=[-44.3281, 6.2470]
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
g10 <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
dat <- raw %>% filter(!if_all(all_of(g10), ~ .x==0))
cat("n_after:", nrow(dat), "\n")
for (c in c("category_one_defects","category_two_defects","quakers")){
x<-dat[[c]]; th<-unname(quantile(x,0.98,na.rm=TRUE)); cat(c,"q98:",th," n_flag(>):",sum(x>th,na.rm=TRUE),"\n")
}' 2>&1n_after: 1338 category_one_defects q98: 5 n_flag(>): 20 category_two_defects q98: 20 n_flag(>): 23 quakers q98: 2 n_flag(>): 25
command: create
path: /app/analysis.R
file_text: #!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-aware rewrite
# ===================================================
# Each column is summarised with a method that matches its distribution instead
# of a single off-the-shelf Tukey/Mahalanobis/rstandard recipe everywhere.
# * altitude : right-skewed -> Tukey fence on log10 scale (k=3) + unit repair
# * defects : zero-inflated -> upper-tail quantile cut-off (98th pct)
# * grades : 7 SCA flavour attrs only -> classical Mahalanobis, chi-sq .999
# * influence : AIC-selected OLS -> leverage>2p/n AND Cook's D>4/n (both)
# * countries : raw mean vs 10% trimmed mean, rank movement
# The script is self-contained: re-running it from an empty /app/outputs/
# reproduces every artefact.
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: the single withdrawn lot with every grade recorded as 0.
# ---------------------------------------------------------------------------
grades10 <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"uniformity", "clean_cup", "sweetness", "cupper_points")
is_sentinel <- rowSums(raw[grades10] == 0, na.rm = TRUE) == length(grades10)
dat <- raw[!is_sentinel, , drop = FALSE]
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L)) # 0-indexed post-drop
grades7 <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "cupper_points")
# ---------------------------------------------------------------------------
# 2. Altitude: log10 Tukey fence (k = 3) + power-of-ten unit repair.
# ---------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
qa <- quantile(log10(alt_pos), c(0.25, 0.75))
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 # back-transform bounds to meters
hi_m <- 10^hi_log
flag_alt <- !is.na(alt) & alt > 0 & (alt < lo_m | alt > hi_m)
n_flag_alt <- sum(flag_alt)
first_numeric_token <- function(s) {
suppressWarnings(as.numeric(str_extract(s, "[0-9]+(\\.[0-9]+)?")))
}
# For each flagged row, test power-of-ten corrections on the first numeric token
# in priority order: /10, then /100, then as-is. Keep the first in-fence value.
alt_corrected <- rep(NA_real_, nrow(dat))
for (i in which(flag_alt)) {
tok <- first_numeric_token(dat$altitude[i])
if (is.na(tok)) next
for (cand in c(tok / 10, tok / 100, tok)) {
if (!is.na(cand) && cand >= lo_m && cand <= hi_m) { alt_corrected[i] <- cand; break }
}
}
n_unit_corrected <- sum(!is.na(alt_corrected))
# altitude used for the regression: recovered meters where available, else raw
alt_for_reg <- ifelse(!is.na(alt_corrected), alt_corrected, alt)
# ---------------------------------------------------------------------------
# 3. Defects: zero-inflated counts -> upper-tail quantile cut-off (98th pct).
# A raw IQR fence collapses (Q1=Q3=0) and flags every non-zero count, so we
# instead threshold the sparse upper tail directly.
# ---------------------------------------------------------------------------
defect_tail <- function(x, p = 0.98) {
upper <- unname(quantile(x, p, na.rm = TRUE))
flagged <- !is.na(x) & x > upper
list(upper_threshold = upper, n = sum(flagged), flag_vec = flagged)
}
d_cat1 <- defect_tail(dat$category_one_defects)
d_cat2 <- defect_tail(dat$category_two_defects)
d_quak <- defect_tail(dat$quakers)
flag_defect <- d_cat1$flag_vec | d_cat2$flag_vec | d_quak$flag_vec
# --------
…[truncated]Existing file was successfully overwritten at /app/analysis.R
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1
Audit complete: n_input=1339, n_after=1338, selected=grades_defects altitude fence [357.45, 4923.82] m, flagged=51, unit-corrected=7 mv flagged=40, influence both=35, composite=105 delta_pct=-20.748 CI=[-44.328, 6.247] (B_ok=500)
cd /app && cat outputs/outlier_report.json
{
"n_input": 1339,
"n_after_zero_score_drop": 1338,
"methods": {
"altitude": "Tukey/IQR fence (k=3) on log10(altitude_mean_meters>0), bounds back-transformed to meters; flagged rows repaired by first in-fence power-of-ten correction (/10,/100,as-is) of the first numeric token in raw altitude",
"defects": "Per-column upper-tail cut-off at the 98th percentile of each zero-inflated count (raw IQR fence collapses since Q1=Q3=0); flag count strictly above threshold, union across columns",
"multivariate_grades": "Classical Mahalanobis distance (sample mean/covariance) on the 7 primary SCA flavour attributes only; chi-squared 99.9% upper-tail threshold with df=7",
"influence": "Lowest finite-AIC OLS among grades_only/grades_altitude/grades_defects (common n); flag leverage>2p/n AND Cook's D>4/n on the selected fit",
"country_ranking": "Countries with n>=10 lots: raw mean vs 10% trimmed mean of total_cup_points, compared by descending-rank movement"
},
"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,
"tr
…[truncated]cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(jsonlite)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
cat("CSV rows:", nrow(csv), " row_id range:", min(csv$row_id), max(csv$row_id), "\n")
cat("row_id sequential 0..n-1:", all(csv$row_id == 0:(nrow(csv)-1)), "\n")
# reload lm and recompute
fit <- readRDS("outputs/coffee_ols_selected.rds")
cat("fit class:", class(fit), " terms:", paste(names(coef(fit)),collapse=","), "\n")
n<-nobs(fit); p<-length(coef(fit))
lev<-hatvalues(fit); ck<-cooks.distance(fit)
hl <- lev > 2*p/n; hc <- ck > 4/n
# all complete -> fit rows map 1:1 to csv rows
cat("recomputed high_lev matches CSV:", all(hl == csv$high_leverage_flag), "\n")
cat("recomputed high_cook matches CSV:", all(hc == csv$high_cook_flag), "\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:", all(comp == csv$composite_drop_flag), "\n")
# defect independent of composite
cat("defect NOT in composite? rows defect-only excluded:", sum(csv$defect_outlier_flag & !comp), "defect rows not composite\n")
# altitude_corrected_m only non-NA when altitude flagged
cat("corrected rows all flagged:", all(csv$altitude_outlier_flag[!is.na(csv$altitude_corrected_m)]), "\n")
cat("n corrected non-NA:", sum(!is.na(csv$altitude_corrected_m)), "\n")
# cross-check JSON top10 against recomputed cook order
js <- fromJSON("outputs/outlier_report.json")
ord <- order(ck, decreasing=TRUE)[1:10]-1L
cat("top10 cook ids match JSON:", all(ord == js$influence$top10_cook_row_ids), "\n")
' 2>&1CSV rows: 1338 row_id range: 0 1337 row_id sequential 0..n-1: TRUE fit class: lm terms: (Intercept),aroma,flavor,aftertaste,acidity,body,balance,cupper_points,ld1,ld2,ldq recomputed high_lev matches CSV: TRUE recomputed high_cook matches CSV: TRUE composite matches CSV: TRUE defect NOT in composite? rows defect-only excluded: 47 defect rows not composite corrected rows all flagged: TRUE n corrected non-NA: 7 top10 cook ids match JSON: TRUE
cd /app && Rscript -e '
suppressPackageStartupMessages({library(jsonlite)})
js <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
errs <- character(0)
need <- function(path, cond) if(!isTRUE(cond)) errs<<-c(errs, path)
# top-level
for (k in c("n_input","n_after_zero_score_drop","methods","altitude","defects","multivariate_grades","influence","country_ranking","sensitivity","flag_interactions")) need(k, k %in% names(js))
for (k in c("altitude","defects","multivariate_grades","influence","country_ranking")) need(paste0("methods.",k), is.character(js$methods[[k]]))
for (k in c("lower_bound_m","upper_bound_m","n_flagged","n_unit_corrected")) need(paste0("altitude.",k), is.numeric(js$altitude[[k]]))
for (c in c("category_one_defects","category_two_defects","quakers")) for(k in c("upper_threshold","n_flagged")) need(paste0("defects.",c,".",k), is.numeric(js$defects[[c]][[k]]))
need("mv.columns_used", length(js$multivariate_grades$columns_used)==7)
for (k in c("threshold","n_flagged")) need(paste0("mv.",k), is.numeric(js$multivariate_grades[[k]]))
for (k in 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")) need(paste0("influence.",k), k %in% names(js$influence))
need("influence.selected_model in set", js$influence$selected_model %in% c("grades_only","grades_altitude","grades_defects"))
need("influence.top10 length10", length(js$influence$top10_cook_row_ids)==10)
for (k in c("n_countries_evaluated","n_rank_changes_ge_2","max_abs_rank_change","per_country")) need(paste0("country.",k), k %in% names(js$country_ranking))
pc <- js$country_ranking$per_country[[1]]
for (k in c("country","n","raw_mean","trimmed_mean","raw_rank","trimmed_rank")) need(paste0("per_country.",k), k %in% names(pc))
for (k in c("coef_flavor_full","coef_flavor_dropped","delta_pct","delta_pct_ci_lower","delta_pct_ci_upper","n_bootstrap_replicates")) need(paste0("sensitivity.",k), is.numeric(js$sensitivity[[k]]))
for (k in 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")) need(paste0("flag_interactions.",k), is.numeric(js$flag_interactions[[k]]))
if(length(errs)==0) cat("SCHEMA OK: all required keys present with expected types\n") else {cat("SCHEMA ERRORS:\n"); print(errs)}
cat("per_country entries:", length(js$country_ranking$per_country), "\n")
' 2>&1SCHEMA OK: all required keys present with expected types per_country entries: 21
cd /app && head -1 outputs/outlier_flags.csv && echo "---expected---" && echo "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" && echo "---corrected rows---" && Rscript -e '
suppressPackageStartupMessages(library(readr))
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
print(as.data.frame(csv[!is.na(csv$altitude_corrected_m), c("row_id","altitude_outlier_flag","altitude_corrected_m","composite_drop_flag")]))
cat("flag counts: alt=",sum(csv$altitude_outlier_flag)," defect=",sum(csv$defect_outlier_flag)," mv=",sum(csv$multivariate_grade_outlier_flag)," lev=",sum(csv$high_leverage_flag)," cook=",sum(csv$high_cook_flag)," composite=",sum(csv$composite_drop_flag),"\n")
' 2>&1row_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 ---expected--- 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 ---corrected rows--- row_id altitude_outlier_flag altitude_corrected_m composite_drop_flag 1 543 TRUE 1100.00 TRUE 2 584 TRUE 518.00 TRUE 3 677 TRUE 518.00 TRUE 4 833 TRUE 518.00 TRUE 5 896 TRUE 1901.64 TRUE 6 1040 TRUE 1100.00 TRUE 7 1144 TRUE 1901.64 TRUE flag counts: alt= 51 defect= 63 mv= 40 lev= 110 cook= 66 composite= 105
command: create path: /app/outputs/audit_memo.md file_text: # Coffee Quality Outlier Audit , Method 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 a clean `/app/outputs/`. ## Data prep The draft dropped lots with `total_cup_points > 0`, which happens to remove the right row but for the wrong reason. The documented sentinel is a *withdrawn submission with every grade recorded as zero*. I drop on that explicit condition (all ten grade columns equal zero), which is self-documenting and robust to any future lot that scores a literal zero total without being blank. Input **1339 → 1338** after the drop. Row ids are 0-indexed in post-drop input order so every artifact aligns. ## Altitude `altitude_mean_meters` is strongly right-skewed (median ~1310 m, max 190,164 m). A Tukey fence on the raw meter scale, as in the draft, is dragged by that tail and produces asymmetric, misleading bounds. I build the fence on `log10(altitude)` over positive values with `k = 3`, then back-transform: **[357.45 m, 4923.82 m]**. 51 rows fall outside. Many are decimal-displacement typos in the raw `altitude` string, so for each flagged row I test power-of-ten corrections on the first numeric token in priority order `/10 → /100 → as-is` and keep the first candidate inside the fence. That repairs **7** rows (e.g. `190164`→1901.64, `11000 metros`→1100); the rest stay `NA`. The draft never attempted recovery at all. ## Defects `category_one_defects`, `category_two_defects` and `quakers` are mass-at-zero counts (85%, 28%, 93% zeros). Their quartiles are `Q1 = Q3 = 0`, so the draft's IQR fence has width zero and flags *every* non-zero count as an outlier , useless for isolating genuinely extreme lots. Instead I cut each column at its **98th percentile** and flag counts strictly above it, directly targeting the sparse upper tail the task describes. Thresholds 5 / 20 / 2 flag 20 / 23 / 25 lots (~1.5–1.9% each); a lot is a defect outlier if any column trips (63 lots). ## Multivariate grades The draft ran Mahalanobis on all ten grade columns. Three of those (`uniformity`, `clean_cup`, `sweetness`) sit at a near-constant 10, so the covariance matrix is near-singular and needed a `tol = 1e-30` hack; distances are then dominated by trivial deviations in degenerate directions. I restrict to the **seven primary SCA flavor attributes** (aroma, flavor, aftertaste, acidity, body, balance, cupper_points), compute a classical Mahalanobis distance from the sample mean and covariance, and threshold at the chi-squared 99.9% upper tail with `df = 7` (**24.32**). This flags **40** joint outliers on a well-conditioned covariance. ## Influence `abs(rstandard) > 2` measures residual size, not influence, and flags ~5% of rows by construction regardless of leverage. I instead fit three OLS candidates on the rows complete in the seven grades and select the lowest **finite** AIC (imputation keeps `n = 1338` identical, so AIC is comparable): grades_only 4459.36, grades_altitude 4461.29, **grades_defects 4418.68** , the defect spec wins. On that fit I flag high leverage (`> 2p/n`, threshold 0.0164) and high Cook's D (`> 4/n`, threshold 0.00299) **separately**, and define influence as **both** (leverage *and* Cook's D): 110 high-leverage, 66 high-Cook, **35** both. The saved `lm` lets the verifier recompute both exactly. ## Country ranking Raw means alone (the draft) are sensitive to one anomalous lot in a small country. For countries with **≥10 lots** (21 of them) I compare the raw mean of `total_cup_points` to a **10% trimmed mean** and rank both descending. **2** countries move ≥2 positions (max change 2): Taiwan slips 14→16 and Honduras rises 20→18 once extreme lots are down-weighted , exactly the fragility a raw mean hides. The per-country breakdown with both means and ranks is in the JSON. ## Sensitivity …[truncated]
New file created successfully at /app/outputs/audit_memo.md
cd /app && echo "word count:" && wc -w outputs/audit_memo.md && echo "--- verify cited interaction numbers ---" && Rscript -e '
suppressPackageStartupMessages(library(jsonlite))
js <- fromJSON("outputs/outlier_report.json")
fi <- js$flag_interactions
cat("0/1/2/3/4:",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,"\n")
cat("alt_only/defect_only/mv_only/infl_only:",fi$n_alt_only,fi$n_defect_only,fi$n_mv_only,fi$n_influence_only,"\n")
cat("mv_and_influence:",fi$n_mv_and_influence,"\n")
cat("sum check:",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,"\n")
'word count: 818 outputs/audit_memo.md --- verify cited interaction numbers --- 0/1/2/3/4: 1186 122 23 7 0 alt_only/defect_only/mv_only/infl_only: 42 47 23 10 mv_and_influence: 13 sum check: 1338
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 = NA)
new_str: write_json(results, file.path(OUTPUT_DIR, "outlier_report.json"),
auto_unbox = FALSE, pretty = TRUE, na = "null", digits = NA)
# ---------------------------------------------------------------------------
# Audit memo. Generated here (not a static file) so the script alone rebuilds
# every artifact on a clean /app/outputs/. Key figures are interpolated from the
# computed objects above so the prose can never drift from the JSON/CSV.
# ---------------------------------------------------------------------------
z1 <- 100 * mean(dat$category_one_defects == 0, na.rm = TRUE)
z2 <- 100 * mean(dat$category_two_defects == 0, na.rm = TRUE)
zq <- 100 * mean(dat$quakers == 0, na.rm = TRUE)
pct_flag_d1 <- 100 * d_cat1$n / n_after
pct_flag_d2 <- 100 * d_cat2$n / n_after
pct_flag_dq <- 100 * d_quak$n / n_after
n_defect_union <- sum(flag_defect)
movers <- country_tbl %>% filter(abs_rank_change >= 2)
movers_txt <- if (nrow(movers) > 0)
paste(sprintf("%s (%d\u2192%d)", movers$country, movers$raw_rank, movers$trimmed_rank),
collapse = ", ") else "none"
memo <- sprintf('# Coffee Quality Outlier Audit \u2014 Method 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 a clean `/app/outputs/`.
## Data prep
The draft dropped lots with `total_cup_points > 0`, which happens to remove the
right row but for the wrong reason. The documented sentinel is a *withdrawn
submission with every grade recorded as zero*. I drop on that explicit
condition (all ten grade columns equal zero), which is self-documenting and
robust to any future lot that scores a literal zero total without being blank.
Input **%d \u2192 %d** after the drop. Row ids are 0-indexed in post-drop input
order so every artifact aligns.
## Altitude
`altitude_mean_meters` is strongly right-skewed (median ~1310 m, max 190,164 m).
A Tukey fence on the raw meter scale, as in the draft, is dragged by that tail
and produces misleading bounds. I build the fence on `log10(altitude)` over
positive values with `k = 3`, then back-transform: **[%.2f m, %.2f m]**. %d rows
fall outside. Many are decimal-displacement typos in the raw `altitude` string,
so for each flagged row I test power-of-ten corrections on the first numeric
token in priority order `/10 \u2192 /100 \u2192 as-is` and keep the first candidate
inside the fence. That repairs **%d** rows (e.g. `190164`\u21921901.64,
`11000 metros`\u21921100); the rest stay `NA`. The draft never attempted recovery.
## Defects
`category_one_defects`, `category_two_defects` and `quakers` are mass-at-zero
counts (%.0f%%, %.0f%%, %.0f%% zeros). Their quartiles are `Q1 = Q3 = 0`, so the
draft'"'"'s IQR fence has width zero and flags *every* non-zero count as an
outlier \u2014 useless for isolating genuinely extreme lots. Instead I cut each
column at its **98th percentile** and flag counts strictly above it, directly
targeting the sparse upper tail the task describes. Thresholds %g / %g / %g flag
%d / %d / %d lots (~%.1f\u2013%.1f%% each); a lot is a defect outlier if any column
trips (%d lots).
## Multivariate grades
The draft ran Mahalanobis on all ten grade columns. Three of those
(`uniformity`, `clean_cup`, `sweetness`) sit at a near-constant 10, so the
covariance is near-singular (it needed a `tol = 1e-30` hack) and distances are
dominated by trivial deviations in degenerate directions. I restrict to the
**seven primary SCA flavor attributes** (aroma, flavor, aftertaste, acidity,
body, balance, cupper_points), compute a classical Mahalanobis distance from the
sample mean and covariance, and threshold at the chi-squared 99.9%% upper tail
with `df = 7` (**%.2f**). This flag
…[truncated]The file /app/analysis.R has been edited successfully.
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1 && echo "=== FILES ===" && ls -la outputs/ && echo "=== MEMO WORDS ===" && wc -w outputs/audit_memo.md && echo "=== HEADINGS ===" && grep '^## ' outputs/audit_memo.md
Error: unexpected string constant in: "The draft applied one off-the-shelf recipe to every column. Each block below states why that recipe fails for the column'"'"" Execution halted
cd /app && sed -i "s/'\"'\"'/\\\\'/g" analysis.R && echo "=== check replacements ===" && grep -n "column.s actual\|draft.s IQR\|flavor.s weight\|Cook.s D" analysis.R
=== check replacements ===
9:# * influence : AIC-selected OLS -> leverage>2p/n AND Cook's D>4/n (both)
161:# top-10 largest Cook's D -> 0-indexed post-drop row ids
287: influence = unbox("Lowest finite-AIC OLS among grades_only/grades_altitude/grades_defects (common n); flag leverage>2p/n AND Cook's D>4/n on the selected fit"),cd /app && grep -n "column\\\\'s\|draft\\\\'s\|flavor\\\\'s\|Cook\\\\'s D\`" analysis.R | head; echo "--- show the escaped apostrophes in memo block ---"; grep -n "\\\\'s" analysis.R
369:states why that recipe fails for the column\'s actual shape and what replaced it. 398:draft\'s IQR fence has width zero and flags *every* non-zero count as an 444:materially inflate flavor\'s weight. A percentile bootstrap --- show the escaped apostrophes in memo block --- 369:states why that recipe fails for the column\'s actual shape and what replaced it. 398:draft\'s IQR fence has width zero and flags *every* non-zero count as an 424:fit I flag high leverage (`> 2p/n`, threshold %.4f) and high Cook\'s D 444:materially inflate flavor\'s weight. A percentile bootstrap
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1 && echo "=== FILES ===" && ls -la outputs/ && echo "=== MEMO WORDS ===" && wc -w outputs/audit_memo.md && echo "=== HEADINGS ===" && grep '^## ' outputs/audit_memo.md
Audit complete: n_input=1339, n_after=1338, selected=grades_defects altitude fence [357.45, 4923.82] m, flagged=51, unit-corrected=7 mv flagged=40, influence both=35, composite=105 delta_pct=-20.748 CI=[-44.328, 6.247] (B_ok=500) === FILES === total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:24 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:24 .. -rw-r--r-- 1 h2tagent h2tagent 5298 Jun 20 18:24 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137127 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 6830 Jun 20 18:24 outlier_report.json === MEMO WORDS === 804 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 , Method 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 a clean `/app/outputs/`. ## Data prep The draft dropped lots with `total_cup_points > 0`, which happens to remove the right row but for the wrong reason. The documented sentinel is a *withdrawn submission with every grade recorded as zero*. I drop on that explicit condition (all ten grade columns equal zero), which is self-documenting and robust to any future lot that scores a literal zero total without being blank. Input **1339 → 1338** after the drop. Row ids are 0-indexed in post-drop input order so every artifact aligns. ## Altitude `altitude_mean_meters` is strongly right-skewed (median ~1310 m, max 190,164 m). A Tukey fence on the raw meter scale, as in the draft, is dragged by that tail and produces misleading bounds. I build the fence on `log10(altitude)` over positive values with `k = 3`, then back-transform: **[357.45 m, 4923.82 m]**. 51 rows fall outside. Many are decimal-displacement typos in the raw `altitude` string, so for each flagged row I test power-of-ten corrections on the first numeric token in priority order `/10 → /100 → as-is` and keep the first candidate inside the fence. That repairs **7** rows (e.g. `190164`→1901.64, `11000 metros`→1100); the rest stay `NA`. The draft never attempted recovery. ## Defects `category_one_defects`, `category_two_defects` and `quakers` are mass-at-zero counts (85%, 28%, 93% zeros). Their quartiles are `Q1 = Q3 = 0`, so the draft's IQR fence has width zero and flags *every* non-zero count as an outlier , useless for isolating genuinely extreme lots. Instead I cut each column at its **98th percentile** and flag counts strictly above it, directly targeting the sparse upper tail the task describes. Thresholds 5 / 20 / 2 flag 20 / 23 / 25 lots (~1.5–1.9% each); a lot is a defect outlier if any column trips (63 lots). ## Multivariate grades The draft ran Mahalanobis on all ten grade columns. Three of those (`uniformity`, `clean_cup`, `sweetness`) sit at a near-constant 10, so the covariance is near-singular (it needed a `tol = 1e-30` hack) and distances are dominated by trivial deviations in degenerate directions. I restrict to the **seven primary SCA flavor attributes** (aroma, flavor, aftertaste, acidity, body, balance, cupper_points), compute a classical Mahalanobis distance from the sample mean and covariance, and threshold at the chi-squared 99.9% upper tail with `df = 7` (**24.32**). This flags **40** joint outliers on a well-conditioned covariance. ## Influence `abs(rstandard) > 2` measures residual size, not influence, and flags ~5% of rows by construction regardless of leverage. I instead fit three OLS candidates on the rows complete in the seven grades and select the lowest **finite** AIC (imputation keeps `n = 1338` identical, so AIC is comparable): grades_only 4459.36, grades_altitude 4461.29, **grades_defects 4418.68** , the defect spec wins. On that fit I flag high leverage (`> 2p/n`, threshold 0.0164) and high Cook's D (`> 4/n`, threshold 0.00299) **separately**, and define influence as **both** (leverage *and* Cook): 110 high-leverage, 66 high-Cook, **35** both. The saved `lm` lets the verifier recompute both exactly. ## Country ranking Raw means alone (the draft) are sensitive to one anomalous lot in a small country. For countries with **≥10 lots** (21 of them) I compare the raw mean of `total_cup_points` to a **10% trimmed mean** and rank both descending. **2** countries move ≥2 positions (max change 2): Taiwan (14→16), Honduras (20→18) , exactly the fragility a raw mean hides. The per-country breakdown with both means and ranks is in the JSON. ## Sensitivity The draft reused an `rstandard`-based composite and never quantified uncertainty. I define composite drop as **altitude OR multivariate …[truncated]
cd /app && cp outputs/outlier_report.json /tmp/run1.json && cp outputs/outlier_flags.csv /tmp/run1.csv && cp outputs/audit_memo.md /tmp/run1.md && 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 "MEMO identical:" && diff -q /tmp/run1.md outputs/audit_memo.md && echo "ALL DETERMINISTIC"
JSON identical: CSV identical: MEMO identical: ALL DETERMINISTIC
cd /tmp && Rscript -e '
# Fresh session, no access to /app data objects -- mimic verifier reloading RDS
fit <- readRDS("/app/outputs/coffee_ols_selected.rds")
n <- nobs(fit); p <- length(coef(fit))
lev <- hatvalues(fit); ck <- cooks.distance(fit)
cat("Reloaded standalone 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("top10 cook (0-idx):",paste(order(ck,decreasing=TRUE)[1:10]-1L,collapse=","),"\n")
cat("formula:",deparse(formula(fit)),"\n")
' 2>&1Reloaded standalone OK. n= 1338 p= 11 n_high_lev: 110 n_high_cook: 66 n_both: 35 top10 cook (0-idx): 1309,1308,1333,1295,1306,1260,1307,1302,1291,1259 formula: total_cup_points ~ aroma + flavor + aftertaste + acidity + body + balance + cupper_points + ld1 + ld2 + ldq
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr);library(dplyr);library(stringr);library(jsonlite)})
js<-fromJSON("outputs/outlier_report.json"); csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
raw<-read_csv("data/coffee_ratings.csv",show_col_types=FALSE)
g10<-c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
dat<-raw[rowSums(raw[g10]==0,na.rm=TRUE)!=10,]
ok<-function(lbl,c) cat(sprintf("[%s] %s\n", ifelse(isTRUE(c),"PASS","FAIL"), lbl))
# counts
ok("n_input=1339", js$n_input==1339 && nrow(raw)==1339)
ok("n_after=1338", js$n_after_zero_score_drop==1338 && nrow(dat)==1338 && nrow(csv)==1338)
# altitude recompute
alt<-dat$altitude_mean_meters; ap<-alt[!is.na(alt)&alt>0]; q<-quantile(log10(ap),c(.25,.75)); iqr<-q[[2]]-q[[1]]
lo<-10^(q[[1]]-3*iqr); hi<-10^(q[[2]]+3*iqr)
fa<-!is.na(alt)&alt>0&(alt<lo|alt>hi)
ok("altitude bounds match JSON", isTRUE(all.equal(lo,js$altitude$lower_bound_m)) && isTRUE(all.equal(hi,js$altitude$upper_bound_m)))
ok("altitude flag matches CSV & JSON n", all(fa==csv$altitude_outlier_flag) && sum(fa)==js$altitude$n_flagged)
# defects recompute
d<-function(x){th<-unname(quantile(x,.98,na.rm=TRUE));list(th=th,f=!is.na(x)&x>th)}
d1<-d(dat$category_one_defects);d2<-d(dat$category_two_defects);dq<-d(dat$quakers)
ok("defect thresholds match JSON", d1$th==js$defects$category_one_defects$upper_threshold & d2$th==js$defects$category_two_defects$upper_threshold & dq$th==js$defects$quakers$upper_threshold)
ok("defect union matches CSV", all((d1$f|d2$f|dq$f)==csv$defect_outlier_flag))
# mahalanobis recompute
g7<-c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G<-as.matrix(dat[,g7]); md<-mahalanobis(G,colMeans(G),cov(G)); fm<-md>qchisq(0.999,7)
ok("mv flag matches CSV & JSON", all(fm==csv$multivariate_grade_outlier_flag) & sum(fm)==js$multivariate_grades$n_flagged & js$multivariate_grades$threshold==qchisq(0.999,7))
# composite
comp<-csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag)
ok("composite matches CSV", all(comp==csv$composite_drop_flag))
ok("defects excluded from composite (some defect-only exist)", sum(csv$defect_outlier_flag & !csv$composite_drop_flag)>0)
# flag interactions recompute from CSV
A<-csv$altitude_outlier_flag;M<-csv$multivariate_grade_outlier_flag;D<-csv$defect_outlier_flag;I<-csv$high_leverage_flag&csv$high_cook_flag
tf<-A+M+D+I; fi<-js$flag_interactions
ok("interaction 0-4 counts", fi$n_with_zero_flags==sum(tf==0)&fi$n_with_one_flag==sum(tf==1)&fi$n_with_two_flags==sum(tf==2)&fi$n_with_three_flags==sum(tf==3)&fi$n_with_four_flags==sum(tf==4))
ok("pairwise intersections", fi$n_alt_and_mv==sum(A&M)&fi$n_alt_and_influence==sum(A&I)&fi$n_mv_and_influence==sum(M&I)&fi$n_alt_mv_influence==sum(A&M&I))
ok("only-counts", fi$n_alt_only==sum(A&!M&!I&!D)&fi$n_mv_only==sum(M&!A&!I&!D)&fi$n_influence_only==sum(I&!A&!M&!D)&fi$n_defect_only==sum(D&!A&!M&!I))
# selected model + AIC ordering
ok("selected grades_defects has min AIC", js$influence$selected_model=="grades_defects" & js$influence$aic_grades_defects<js$influence$aic_grades_only & js$influence$aic_grades_defects<js$influence$aic_grades_altitude)
# country ranks internal consistency
pc<-js$country_ranking$per_country
ok("country raw_rank consistent with raw_mean order", all(rank(-pc$raw_mean,ties.method="min")==pc$raw_rank))
ok("country trimmed_rank consistent", all(rank(-pc$trimmed_mean,ties.method="min")==pc$trimmed_rank))
ok("n_rank_changes_ge_2 matches", sum(abs(pc$raw_rank-pc$trimmed_rank)>=2)==js$country_ranking$n_rank_changes_ge_2)
ok("all countries n>=10", all(pc$n>=10))
# sensitivity delta sign
ok("delta_pct = 100*(drop-full)/full", isTRUE(all.equal(js$sensitivity$delta_pct, 100*(js$sensitivity$coef_flavor_dropped-js$sensitivity$coef_flavor_full)/js$sensitivity$coef_flavor_full)))
ok("CI brackets ordering", js$sensitivity$delta_pct_ci_lower<=js$sensitivity$delta_p
…[truncated][PASS] n_input=1339 [PASS] n_after=1338 [PASS] altitude bounds match JSON [PASS] altitude flag matches CSV & JSON n [PASS] defect thresholds match JSON [PASS] defect union matches CSV [FAIL] mv flag matches CSV & JSON [PASS] composite matches CSV [PASS] defects excluded from composite (some defect-only exist) [PASS] interaction 0-4 counts [PASS] pairwise intersections [PASS] only-counts [PASS] selected grades_defects has min AIC [PASS] country raw_rank consistent with raw_mean order [PASS] country trimmed_rank consistent [PASS] n_rank_changes_ge_2 matches [PASS] all countries n>=10 [PASS] delta_pct = 100*(drop-full)/full [PASS] CI brackets ordering
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr);library(dplyr);library(jsonlite)})
js<-fromJSON("outputs/outlier_report.json"); csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
raw<-read_csv("data/coffee_ratings.csv",show_col_types=FALSE)
g10<-c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
dat<-raw[rowSums(raw[g10]==0,na.rm=TRUE)!=10,]
g7<-c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G<-as.matrix(dat[,g7]); md<-mahalanobis(G,colMeans(G),cov(G)); fm<-md>qchisq(0.999,7)
cat("flag vec matches CSV:", all(fm==csv$multivariate_grade_outlier_flag), "\n")
cat("sum fm:", sum(fm), " JSON n_flagged:", js$multivariate_grades$n_flagged, "\n")
cat("threshold JSON:", sprintf("%.15f",js$multivariate_grades$threshold), " recompute:", sprintf("%.15f",qchisq(0.999,7)), "\n")
cat("threshold equal (==):", js$multivariate_grades$threshold==qchisq(0.999,7), " all.equal:", isTRUE(all.equal(js$multivariate_grades$threshold,qchisq(0.999,7))), "\n")
' 2>&1flag vec matches CSV: TRUE sum fm: 40 JSON n_flagged: 40 threshold JSON: 24.321886347856900 recompute: 24.321886347856854 threshold equal (==): FALSE all.equal: TRUE
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr);library(dplyr);library(jsonlite)})
js<-fromJSON("outputs/outlier_report.json"); csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
raw<-read_csv("data/coffee_ratings.csv",show_col_types=FALSE)
g10<-c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
dat<-raw[rowSums(raw[g10]==0,na.rm=TRUE)!=10,]
g7<-c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G<-as.matrix(dat[,g7]); md<-mahalanobis(G,colMeans(G),cov(G)); fm<-md>qchisq(0.999,7)
feq<-function(a,b,tol=1e-6) abs(a-b)<tol
ok<-function(lbl,c) cat(sprintf("[%s] %s\n", ifelse(isTRUE(c),"PASS","FAIL"), lbl))
ok("mv flag matches CSV & count & threshold(tol)", all(fm==csv$multivariate_grade_outlier_flag) & sum(fm)==js$multivariate_grades$n_flagged & feq(js$multivariate_grades$threshold,qchisq(0.999,7)))
ok("columns_used == 7 primary attrs in order", identical(js$multivariate_grades$columns_used, g7))
ok("leverage_threshold = 2p/n (tol)", feq(js$influence$leverage_threshold, 2*11/1338))
ok("cook_threshold = 4/n (tol)", feq(js$influence$cook_threshold, 4/1338))
ok("n_bootstrap_replicates==500", js$sensitivity$n_bootstrap_replicates==500)
ok("altitude n_unit_corrected==7", js$altitude$n_unit_corrected==7)
ok("corrected rows are subset of flagged rows", all(csv$altitude_outlier_flag[!is.na(csv$altitude_corrected_m)]))
' 2>&1[PASS] mv flag matches CSV & count & threshold(tol) [PASS] columns_used == 7 primary attrs in order [PASS] leverage_threshold = 2p/n (tol) [PASS] cook_threshold = 4/n (tol) [PASS] n_bootstrap_replicates==500 [PASS] altitude n_unit_corrected==7 [PASS] corrected rows are subset of flagged rows
cd /app && rm -rf outputs && mkdir outputs && Rscript -e 'source("analysis.R")' 2>&1 | head -30; echo "=== exit: $? ==="; echo "--- any R warnings captured ---"; Rscript -e 'options(warn=1); source("/app/analysis.R")' 2>&1 | grep -i "warn" | headAudit complete: n_input=1339, n_after=1338, selected=grades_defects altitude fence [357.45, 4923.82] m, flagged=51, unit-corrected=7 mv flagged=40, influence both=35, composite=105 delta_pct=-20.748 CI=[-44.328, 6.247] (B_ok=500) === exit: 0 === --- any R warnings captured ---
cd /app && sed -n '1,120p' analysis.R
#!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-aware rewrite
# ===================================================
# Each column is summarised with a method that matches its distribution instead
# of a single off-the-shelf Tukey/Mahalanobis/rstandard recipe everywhere.
# * altitude : right-skewed -> Tukey fence on log10 scale (k=3) + unit repair
# * defects : zero-inflated -> upper-tail quantile cut-off (98th pct)
# * grades : 7 SCA flavour attrs only -> classical Mahalanobis, chi-sq .999
# * influence : AIC-selected OLS -> leverage>2p/n AND Cook's D>4/n (both)
# * countries : raw mean vs 10% trimmed mean, rank movement
# The script is self-contained: re-running it from an empty /app/outputs/
# reproduces every artefact.
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: the single withdrawn lot with every grade recorded as 0.
# ---------------------------------------------------------------------------
grades10 <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"uniformity", "clean_cup", "sweetness", "cupper_points")
is_sentinel <- rowSums(raw[grades10] == 0, na.rm = TRUE) == length(grades10)
dat <- raw[!is_sentinel, , drop = FALSE]
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L)) # 0-indexed post-drop
grades7 <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "cupper_points")
# ---------------------------------------------------------------------------
# 2. Altitude: log10 Tukey fence (k = 3) + power-of-ten unit repair.
# ---------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
qa <- quantile(log10(alt_pos), c(0.25, 0.75))
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 # back-transform bounds to meters
hi_m <- 10^hi_log
flag_alt <- !is.na(alt) & alt > 0 & (alt < lo_m | alt > hi_m)
n_flag_alt <- sum(flag_alt)
first_numeric_token <- function(s) {
suppressWarnings(as.numeric(str_extract(s, "[0-9]+(\\.[0-9]+)?")))
}
# For each flagged row, test power-of-ten corrections on the first numeric token
# in priority order: /10, then /100, then as-is. Keep the first in-fence value.
alt_corrected <- rep(NA_real_, nrow(dat))
for (i in which(flag_alt)) {
tok <- first_numeric_token(dat$altitude[i])
if (is.na(tok)) next
for (cand in c(tok / 10, tok / 100, tok)) {
if (!is.na(cand) && cand >= lo_m && cand <= hi_m) { alt_corrected[i] <- cand; break }
}
}
n_unit_corrected <- sum(!is.na(alt_corrected))
# altitude used for the regression: recovered meters where available, else raw
alt_for_reg <- ifelse(!is.na(alt_corrected), alt_corrected, alt)
# ---------------------------------------------------------------------------
# 3. Defects: zero-inflated counts -> upper-tail quantile cut-off (98th pct).
# A raw IQR fence collapses (Q1=Q3=0) and flags every non-zero count, so we
# instead threshold the sparse upper tail directly.
# ---------------------------------------------------------------------------
defect_tail <- function(x, p = 0.98) {
upper <- unname(quantile(x, p, na.rm = TRUE))
flagged <- !is.na(x) & x > upper
list(upper_threshold = upper, n = sum(flagged), flag_vec = flagged)
}
d_cat1 <- defect_tail(dat$category_one_defects)
d_cat2 <- defect_tail(dat$category_two_defects)
d_quak <- defect_tail(dat$quakers)
flag_defect <- d_cat1$flag_vec | d_cat2$flag_vec | d_quak$flag_vec
# ---------------------------------------------------------
…[truncated]cd /app && sed -n '120,210p' analysis.R
fit_df_alt <- fit_df; fit_df_alt$log_elev <- log_elev
fit_altitude <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points + log_elev, data = fit_df_alt)
# grades_defects: + log1p of the three defect columns, mean-impute rare NA
fit_df_def <- fit_df
fit_df_def$ld1 <- log1p(mean_impute(dat$category_one_defects[G_complete]))
fit_df_def$ld2 <- log1p(mean_impute(dat$category_two_defects[G_complete]))
fit_df_def$ldq <- log1p(mean_impute(dat$quakers[G_complete]))
fit_defects <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points + ld1 + ld2 + ldq,
data = fit_df_def)
aic_vals <- c(grades_only = AIC(fit_grades),
grades_altitude = AIC(fit_altitude),
grades_defects = AIC(fit_defects))
fit_list <- list(grades_only = fit_grades, grades_altitude = fit_altitude,
grades_defects = fit_defects)
finite_aic <- aic_vals[is.finite(aic_vals)]
if (length(finite_aic) == 0) {
selected_model <- "grades_only"
} else {
selected_model <- names(finite_aic)[which.min(finite_aic)]
}
fit_sel <- fit_list[[selected_model]]
n_fit <- nobs(fit_sel)
p_fit <- length(coef(fit_sel))
lev_th <- 2 * p_fit / n_fit
cook_th <- 4 / n_fit
lev_vals <- hatvalues(fit_sel)
cook_vals <- cooks.distance(fit_sel)
high_lev <- rep(FALSE, nrow(dat)); high_lev[G_complete] <- lev_vals > lev_th
high_cook <- rep(FALSE, nrow(dat)); high_cook[G_complete] <- cook_vals > cook_th
infl_both <- high_lev & high_cook
n_high_lev <- sum(high_lev)
n_high_cook <- sum(high_cook)
n_both <- sum(infl_both)
# top-10 largest Cook's D -> 0-indexed post-drop row ids
cc_idx <- which(G_complete)
ord_cook <- order(cook_vals, decreasing = TRUE)[seq_len(min(10, length(cook_vals)))]
top10_row_ids <- as.integer(dat$row_id[cc_idx[ord_cook]])
# ---------------------------------------------------------------------------
# 6. Country ranking: raw mean vs 10% trimmed mean (countries with >= 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.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)
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 of the grades_only `flavor` coefficient.
# Composite = altitude OR multivariate OR (leverage AND Cook). Defects excluded.
# ---------------------------------------------------------------------------
composite <- flag_alt | flag_mv | infl_both
fit_full <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points, data = dat[G_complete, ])
keep_drop <- G_complete & !composite
fit_drop <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points, data = dat[keep_drop, ])
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 delta_pct. Resample post-drop lots; carry the
# ORIGINAL composite flags across replicates (do not recompute per replicate).
boot_df <- dat[G_complete, , drop = FALSE]
comp_boot <- composite[G_complete]
n_boot_rows <- nrow(boot_df)
set.seed(20260512)
B <- 500
grades_only_fit <- function(d) {
tryCatch
…[truncated]cd /app && ls -la outputs/ && echo "=== JSON valid parse ===" && python3 -c "import json; d=json.load(open('outputs/outlier_report.json')); print('keys:', list(d.keys())); print('selected_model:', d['influence']['selected_model']); print('per_country entries:', len(d['country_ranking']['per_country']))" 2>&1total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:25 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:25 .. -rw-r--r-- 1 h2tagent h2tagent 5298 Jun 20 18:25 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137127 Jun 20 18:25 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78269 Jun 20 18:25 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6830 Jun 20 18:25 outlier_report.json === JSON valid parse === keys: ['n_input', 'n_after_zero_score_drop', 'methods', 'altitude', 'defects', 'multivariate_grades', 'influence', 'country_ranking', 'sensitivity', 'flag_interactions'] selected_model: grades_defects per_country entries: 21
cd /app && python3 -c "
import pandas as pd
df = pd.read_csv('outputs/outlier_flags.csv')
print('shape:', df.shape)
print('cols:', list(df.columns))
print('dtypes of flags:', df['altitude_outlier_flag'].dtype, df['composite_drop_flag'].dtype)
print('flag sums:', {c: int(df[c].sum()) for c in ['altitude_outlier_flag','defect_outlier_flag','multivariate_grade_outlier_flag','high_leverage_flag','high_cook_flag','composite_drop_flag']})
print('corrected non-null:', int(df['altitude_corrected_m'].notna().sum()))
print('row_id 0..n-1:', df['row_id'].tolist()==list(range(len(df))))
# composite consistency in python
comp = df['altitude_outlier_flag'] | df['multivariate_grade_outlier_flag'] | (df['high_leverage_flag'] & df['high_cook_flag'])
print('composite matches:', bool((comp==df['composite_drop_flag']).all()))
"shape: (1338, 10)
cols: ['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']
dtypes of flags: bool bool
flag sums: {'altitude_outlier_flag': 51, 'defect_outlier_flag': 63, 'multivariate_grade_outlier_flag': 40, 'high_leverage_flag': 110, 'high_cook_flag': 66, 'composite_drop_flag': 105}
corrected non-null: 7
row_id 0..n-1: True
composite matches: TrueNo file changes, the agent left the workspace untouched.
# SyncVals normalized verifier summary
# task: coffee-ratings-outliers
# attempt: 3 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).
----------------------------------------------------------------------
shape: (1338, 10)
cols: ['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']
dtypes of flags: bool bool
flag sums: {'altitude_outlier_flag': 51, 'defect_outlier_flag': 63, 'multivariate_grade_outlier_flag': 40, 'high_leverage_flag': 110, 'high_cook_flag': 66, 'composite_drop_flag': 105}
corrected non-null: 7
row_id 0..n-1: True
composite matches: TrueReproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_7e4c5961915240f8. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_7e4c5961915240f8 · verifier authoritative; classifier explanatory.