tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.cd /app && ls -la && echo "---DATA---" && ls -la data/ && echo "---OUTPUTS---" && ls -la outputs/ 2>/dev/null || echo "no outputs dir"
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 ---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 ---OUTPUTS--- total 8 drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:46 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 17:46 ..
cd /app && cat analysis.R
#!/usr/bin/env Rscript
# Coffee outlier audit -- first draft
# ===================================
# Quick first pass before the deadline. Loaded the CSV, threw mean / sd / IQR
# at each column, used standardised residuals to flag influential lots, and
# took simple country means. Most of this should be good enough but the
# numbers feel off in a few places.
suppressPackageStartupMessages({
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(jsonlite)
library(purrr)
})
DATA_PATH <- "/app/data/coffee_ratings.csv"
OUTPUT_DIR <- "/app/outputs"
dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE)
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
dat <- raw %>% filter(total_cup_points > 0)
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L))
# ---- Altitude
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
q <- quantile(alt_pos, c(0.25, 0.75))
iqr_alt <- q[[2]] - q[[1]]
lo_m <- q[[1]] - 3 * iqr_alt
hi_m <- q[[2]] + 3 * iqr_alt
flag_alt <- !is.na(alt) & (alt < lo_m | alt > hi_m)
n_flag_alt <- sum(flag_alt)
alt_corrected <- rep(NA_real_, nrow(dat))
unit_corrected <- rep(FALSE, nrow(dat))
n_unit_corrected <- 0L
# ---- Defect counts
defect_summary <- function(x) {
q <- quantile(x, c(0.25, 0.75), na.rm = TRUE)
upper <- q[[2]] + 3 * (q[[2]] - q[[1]])
flagged <- !is.na(x) & x > upper
list(upper_threshold = upper, n = sum(flagged), flag_vec = flagged)
}
d_cat1 <- defect_summary(dat$category_one_defects)
d_cat2 <- defect_summary(dat$category_two_defects)
d_quak <- defect_summary(dat$quakers)
flag_defect <- d_cat1$flag_vec | d_cat2$flag_vec | d_quak$flag_vec
# ---- Multivariate grades (every grade column)
g_all <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "uniformity", "clean_cup",
"sweetness", "cupper_points")
G <- as.matrix(dat[, g_all])
G_complete <- complete.cases(G)
mu <- colMeans(G[G_complete, ])
S <- cov(G[G_complete, ])
md2 <- rep(NA_real_, nrow(G))
md2[G_complete] <- mahalanobis(G[G_complete, ], mu, S, tol = 1e-30)
md2_thresh <- qchisq(0.999, df = length(g_all))
flag_mv <- !is.na(md2) & md2 > md2_thresh
n_flag_mv <- sum(flag_mv)
# ---- Influence: standardised residuals
fit_df <- dat[G_complete, ]
fit <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points, data = fit_df)
rstd <- rstandard(fit)
infl_local <- abs(rstd) > 2
n_fit <- length(rstd)
p_fit <- length(coef(fit))
infl_full <- rep(FALSE, nrow(dat))
infl_full[G_complete] <- infl_local
n_high_lev <- sum(infl_full)
n_high_cook <- sum(infl_full)
n_both <- sum(infl_full)
top_idx <- order(-abs(rstd))[1:10]
top10_row_ids <- as.integer(dat$row_id[which(G_complete)[top_idx]])
# ---- Country ranking (raw means)
country_tbl <- dat %>%
filter(!is.na(country_of_origin)) %>%
group_by(country = country_of_origin) %>%
summarise(n = n(),
raw_mean = mean(total_cup_points),
.groups = "drop") %>%
filter(n >= 10) %>%
arrange(desc(raw_mean)) %>%
mutate(raw_rank = rank(-raw_mean, ties.method = "min"),
trimmed_mean = raw_mean,
trimmed_rank = raw_rank,
abs_rank_change = 0L)
n_countries <- nrow(country_tbl)
n_changes_ge_2 <- 0L
max_change <- 0L
# ---- Composite + sensitivity
composite <- flag_alt | flag_mv | infl_full
keep_for_sens <- !composite & G_complete
fit_drop <- lm(total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points, data = dat[keep_for_sens, ])
coef_full <- unname(coef(fit)["flavor"])
coef_dropped <- unname(coef(fit_drop)["flavor"])
delta_pct <- 100 * (coef_dropped - coef_full) / coef_full
# ---- Flag interactions
A <- flag_alt; M <- flag_mv; D <- flag_defect; I <- infl_full
total_flags <- as.integer(A) + as.integer(M) + as.integer(D) + as.integer(I)
flag_interactions <- list(
n_with_zero_flags = unbox(as.integer(sum(total_flags == 0L))),
…[truncated]cd /app && head -1 data/coffee_ratings.csv | tr ',' '\n' | nl
1 total_cup_points
2 species
3 owner
4 country_of_origin
5 farm_name
6 lot_number
7 mill
8 ico_number
9 company
10 altitude
11 region
12 producer
13 number_of_bags
14 bag_weight
15 in_country_partner
16 harvest_year
17 grading_date
18 owner_1
19 variety
20 processing_method
21 aroma
22 flavor
23 aftertaste
24 acidity
25 body
26 balance
27 uniformity
28 clean_cup
29 sweetness
30 cupper_points
31 moisture
32 category_one_defects
33 quakers
34 color
35 category_two_defects
36 expiration
37 certification_body
38 certification_address
39 certification_contact
40 unit_of_measurement
41 altitude_low_meters
42 altitude_high_meters
43 altitude_mean_meterscd /app && R --version | head -1 && echo "---PKGS---" && Rscript -e 'for (p in c("dplyr","tidyr","readr","stringr","jsonlite","purrr","MASS")) 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 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")
cat("total_cup_points summary:\n"); print(summary(raw$total_cup_points))
cat("\nRows with total_cup_points == 0:\n")
g <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
z <- raw %>% filter(total_cup_points == 0)
cat("count:", nrow(z), "\n")
print(as.data.frame(z[, c("total_cup_points", g)]))
cat("\nRows with any grade==0 (of the seven primary):\n")
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
anyzero <- raw %>% filter(if_any(all_of(g7), ~ .x == 0))
cat("count any primary grade zero:", nrow(anyzero), "\n")
print(as.data.frame(anyzero[, c("total_cup_points", g)]))
'n_input: 1339 total_cup_points summary: Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 81.08 82.50 82.09 83.67 90.58 Rows with total_cup_points == 0: count: 1 total_cup_points aroma flavor aftertaste acidity body balance uniformity 1 0 0 0 0 0 0 0 0 clean_cup sweetness cupper_points 1 0 0 0 Rows with any grade==0 (of the seven primary): count any primary grade zero: 1 total_cup_points aroma flavor aftertaste acidity body balance uniformity 1 0 0 0 0 0 0 0 0 clean_cup sweetness cupper_points 1 0 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)
alt <- dat$altitude_mean_meters
cat("altitude_mean_meters summary:\n"); print(summary(alt))
cat("n NA:", sum(is.na(alt)), " n<=0:", sum(alt<=0, na.rm=TRUE), "\n")
cat("n positive:", sum(alt>0, na.rm=TRUE), "\n\n")
apos <- alt[!is.na(alt) & alt>0]
cat("quantiles of positive altitude:\n"); print(quantile(apos, c(0,.01,.05,.25,.5,.75,.95,.99,1)))
cat("\nlog10 scale quantiles:\n"); print(quantile(log10(apos), c(0,.01,.05,.25,.5,.75,.95,.99,1)))
cat("\nLargest 20 altitudes:\n"); print(sort(apos, decreasing=TRUE)[1:20])
cat("\nSmallest 20 positive altitudes:\n"); print(sort(apos)[1:20])
'altitude_mean_meters summary:
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 1100 1311 1775 1600 190164 230
n NA: 230 n<=0: 0
n positive: 1108
quantiles of positive altitude:
0% 1% 5% 25% 50% 75% 95% 99%
1.00 1.77 442.00 1100.00 1310.64 1600.00 1893.00 3280.00
100%
190164.00
log10 scale quantiles:
0% 1% 5% 25% 50% 75% 95%
0.00000000 0.07554269 2.64542227 3.04139269 3.11748342 3.20411998 3.27714509
99% 100%
3.51587384 5.27912830
Largest 20 altitudes:
[1] 190164 190164 110000 11000 4287 4001 3850 3845 3825 3800
[11] 3500 3280 3280 3280 3170 3140 3000 3000 2560 2560
Smallest 20 positive altitudes:
[1] 1 1 1 1 1 1 1 1 1 1 1 1 12 12 12 13 13 40 50
[20] 100cd /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
apos <- alt[!is.na(alt) & alt>0]
la <- log10(apos)
q <- quantile(la, c(.25,.75))
iqr <- q[[2]]-q[[1]]
lo <- q[[1]] - 3*iqr
hi <- q[[2]] + 3*iqr
cat(sprintf("log10 Q1=%.5f Q3=%.5f IQR=%.5f\n", q[[1]], q[[2]], iqr))
cat(sprintf("log10 fence: [%.5f, %.5f]\n", lo, hi))
cat(sprintf("meter fence: [%.4f, %.4f]\n", 10^lo, 10^hi))
flag <- !is.na(alt) & alt>0 & (log10(alt) < lo | log10(alt) > hi)
# also NA handling: rows with NA altitude can not be flagged
cat("n flagged (positive only):", sum(flag), "\n")
fl <- dat[flag, c("altitude","altitude_low_meters","altitude_high_meters","altitude_mean_meters","unit_of_measurement")]
print(as.data.frame(fl), row.names=FALSE)
'log10 Q1=3.04139 Q3=3.20412 IQR=0.16273
log10 fence: [2.55321, 3.69230]
meter fence: [357.4463, 4923.8167]
n flagged (positive only): 51
altitude altitude_low_meters altitude_high_meters
1.2 12.0000 12.0000
1.2 12.0000 12.0000
-1 1.0000 1.0000
350 350.0000 350.0000
170 170.0000 170.0000
approx. 1'500 metres a.s.l. 1.0000 500.0000
1 1.0000 1.0000
175 175.0000 175.0000
16.8 168.0000 168.0000
approx.1'500 metres a.s.l. 1.0000 500.0000
1 1.0000 1.0000
350 350.0000 350.0000
11000 metros 11000.0000 11000.0000
518 157.8864 157.8864
1.65 165.0000 165.0000
1 1.0000 1.0000
200 m 200.0000 200.0000
518 157.8864 157.8864
1 1.0000 1.0000
200 m 200.0000 200.0000
1 1.0000 1.0000
1 1.0000 1.0000
1 1.0000 1.0000
1.3 13.0000 13.0000
518 157.8864 157.8864
1 1.0000 1.0000
1 1.0000 1.0000
1 1.0000 1.0000
150 m 150.0000 150.0000
190164 190164.0000 190164.0000
1.3 13.0000 13.0000
160-200m 160.0000 200.0000
1.25 125.0000 125.0000
110公尺 110.0000 110.0000
1 1.0000 1.0000
200 200.0000 200.0000
300 m 300.0000 300.0000
100 100.0000 100.0000
1100.00 mosl 110000.0000 110000.0000
1.85 185.0000 185.0000
50 50.0000 50.0000
280 280.0000 280.0000
200 m 200.0000 200.0000
150 m 150.0000 150.0000
250 m 250.0000 250.0000
250 m 250.0000 250.0000
1901.64 190164.0000 190164.0000
200 200.0000 200.0000
12oo 12.0000 12.0000
~350m 350.0000 350.0000
40 40.0000 40.0000
altitude_mean_meters unit_of_measurement
12.0000 m
12.0000 m
1.0000 m
350.0000 m
1
…[truncated]cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(stringr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
alt <- dat$altitude_mean_meters
apos <- alt[!is.na(alt) & alt>0]
la <- log10(apos); q <- quantile(la, c(.25,.75)); iqr<-q[[2]]-q[[1]]
lo_m <- 10^(q[[1]]-3*iqr); hi_m <- 10^(q[[2]]+3*iqr)
flag <- !is.na(alt) & alt>0 & (alt < lo_m | alt > hi_m)
idx <- which(flag)
first_tok <- function(s){ m <- str_extract(s, "[0-9]+(?:\\.[0-9]+)?"); as.numeric(m) }
corrected <- 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_m & cands<=hi_m]
if(length(hit)>0){ corrected[i]<-hit[1]; ncorr<-ncorr+1L }
}
cat("n flagged:", length(idx), " n corrected:", ncorr, "\n\n")
res <- data.frame(altitude=dat$altitude[idx], alt_mean=alt[idx], tok=sapply(dat$altitude[idx], first_tok), corrected=corrected[idx])
print(res, row.names=FALSE)
' 2>&1n flagged: 51 n corrected: 7
altitude alt_mean 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)), " n:", length(x), "\n")
cat("pct zero:", round(mean(x==0, na.rm=TRUE)*100,2), "%\n")
print(table(x, useNA="ifany"))
cat("quantiles:\n"); print(quantile(x, c(.5,.9,.95,.98,.99,.995,1), na.rm=TRUE))
# raw Tukey
q<-quantile(x,c(.25,.75),na.rm=TRUE); up<-q[[2]]+3*(q[[2]]-q[[1]])
cat("raw Tukey k3 upper:", up, " n>up:", sum(x>up,na.rm=TRUE), "\n")
# Q98/Q99 empirical
cat("\n")
}
'==== category_one_defects ====
NA: 0 n: 1338
pct zero: 84.9 %
x
0 1 2 3 4 5 6 7 8 9 10 11 12 15 20 23
1136 101 38 18 16 9 3 3 1 1 4 1 1 1 1 1
31 63
2 1
quantiles:
50% 90% 95% 98% 99% 99.5% 100%
0.000 1.000 2.000 5.000 7.630 11.315 63.000
raw Tukey k3 upper: 0 n>up: 202
==== category_two_defects ====
NA: 0 n: 1338
pct zero: 27.88 %
x
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
373 207 181 134 118 73 44 41 29 23 21 8 13 9 9 6 6 7 2 3
20 21 22 23 24 26 27 28 29 30 31 32 34 38 40 45 47 55
8 2 1 2 1 3 1 1 2 2 1 1 1 1 1 1 1 1
quantiles:
50% 90% 95% 98% 99% 99.5% 100%
2.000 9.000 13.000 20.000 26.630 31.315 55.000
raw Tukey k3 upper: 16 n>up: 43
==== quakers ====
NA: 1 n: 1338
pct zero: 92.97 %
x
0 1 2 3 4 5 6 7 8 9 11 <NA>
1243 39 30 5 5 5 4 3 1 1 1 1
quantiles:
50% 90% 95% 98% 99% 99.5% 100%
0 0 1 2 5 6 11
raw Tukey k3 upper: 0 n>up: 94cd /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]]; thr <- quantile(x, 0.99, na.rm=TRUE, type=7)
fl <- !is.na(x) & x > thr
cat(sprintf("%-22s Q99=%.3f n_flag=%d (%.2f%%)\n", col, thr, sum(fl), 100*mean(fl)))
}
cat("\n--- Multivariate (7 primary grades) ---\n")
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G <- as.matrix(dat[, g7]); cc <- complete.cases(G)
cat("complete rows:", sum(cc), "\n")
mu <- colMeans(G[cc,]); S <- cov(G[cc,])
cat("variances of 7 primary grades:\n"); print(round(diag(S),4))
md2 <- rep(NA_real_, nrow(G)); md2[cc] <- mahalanobis(G[cc,], mu, S)
thr <- qchisq(0.999, df=length(g7))
cat(sprintf("chisq 99.9%% df=7 threshold=%.5f\n", thr))
cat("n flagged:", sum(!is.na(md2) & md2>thr), "\n")
cat("\nFor contrast, variances of near-constant cleanliness cols:\n")
print(round(sapply(dat[,c("uniformity","clean_cup","sweetness")], var, na.rm=TRUE),4))
'category_one_defects Q99=7.630 n_flag=14 (1.05%)
category_two_defects Q99=26.630 n_flag=14 (1.05%)
quakers Q99=5.000 n_flag=10 (0.75%)
--- Multivariate (7 primary grades) ---
complete rows: 1338
variances of 7 primary grades:
aroma flavor aftertaste acidity body
0.0998 0.1165 0.1227 0.1019 0.0947
balance cupper_points
0.1251 0.1822
chisq 99.9% df=7 threshold=24.32189
n flagged: 40
For contrast, variances of near-constant cleanliness cols:
uniformity clean_cup sweetness
0.2354 0.5116 0.3071cd /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 correction
alt <- dat$altitude_mean_meters
apos <- alt[!is.na(alt) & alt>0]; la<-log10(apos); q<-quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]
lo_m<-10^(q[[1]]-3*iqr); hi_m<-10^(q[[2]]+3*iqr)
flag<-!is.na(alt)&alt>0&(alt<lo_m|alt>hi_m)
ft<-function(s){as.numeric(str_extract(s,"[0-9]+(?:\\.[0-9]+)?"))}
corrected<-rep(NA_real_,nrow(dat))
for(i in which(flag)){tok<-ft(dat$altitude[i]); if(is.na(tok))next; cc<-c(tok/10,tok/100,tok); h<-cc[cc>=lo_m&cc<=hi_m]; if(length(h))corrected[i]<-h[1]}
elev <- ifelse(!is.na(corrected), corrected, alt)
le <- log10(elev); le[!is.finite(le)] <- mean(le[is.finite(le)]);
# defects
imp <- function(x){x[is.na(x)]<-mean(x,na.rm=TRUE); x}
cc <- complete.cases(as.matrix(dat[,g7])) # all TRUE
fr <- dat[cc,]; le_f <- le[cc]
f1 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points, data=fr)
fr$le <- le_f
f2 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+le, data=fr)
fr$d1 <- log1p(imp(fr$category_one_defects)); fr$d2<-log1p(imp(fr$category_two_defects)); fr$dq<-log1p(imp(fr$quakers))
f3 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+d1+d2+dq, data=fr)
cat(sprintf("n1=%d n2=%d n3=%d\n", nobs(f1),nobs(f2),nobs(f3)))
cat(sprintf("AIC grades_only = %.6f\n", AIC(f1)))
cat(sprintf("AIC grades_altitude = %.6f\n", AIC(f2)))
cat(sprintf("AIC grades_defects = %.6f\n", AIC(f3)))
aics<-c(grades_only=AIC(f1),grades_altitude=AIC(f2),grades_defects=AIC(f3))
fin<-aics[is.finite(aics)]; cat("selected:", names(fin)[which.min(fin)], "\n")
'n1=1338 n2=1338 n3=1338 AIC grades_only = 4459.361516 AIC grades_altitude = 4461.292330 AIC grades_defects = 4418.680074 selected: grades_defects
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) %>% mutate(row_id=as.integer(seq_len(n())-1L))
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
imp <- function(x){x[is.na(x)]<-mean(x,na.rm=TRUE); x}
fr <- dat
fr$d1<-log1p(imp(fr$category_one_defects)); fr$d2<-log1p(imp(fr$category_two_defects)); fr$dq<-log1p(imp(fr$quakers))
f3 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+d1+d2+dq, data=fr)
n<-nobs(f3); p<-length(coef(f3))
lev<-hatvalues(f3); ck<-cooks.distance(f3)
lev_thr<-2*p/n; ck_thr<-4/n
cat(sprintf("n=%d p=%d lev_thr=%.6f cook_thr=%.6f\n",n,p,lev_thr,ck_thr))
hl<-lev>lev_thr; hc<-ck>ck_thr
cat(sprintf("n_high_leverage=%d n_high_cook=%d n_both=%d\n", sum(hl),sum(hc),sum(hl&hc)))
ord<-order(-ck)[1:10]
cat("top10 cook row_ids:", paste(dat$row_id[ord],collapse=","), "\n")
cat("top10 cook values:", paste(round(ck[ord],4),collapse=","), "\n")
'n=1338 p=11 lev_thr=0.016442 cook_thr=0.002990 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 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)
tb <- 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"),
abs_rank_change=abs(raw_rank-trimmed_rank)) %>%
arrange(raw_rank)
cat("n countries:", nrow(tb), "\n")
cat("n rank changes >=2:", sum(tb$abs_rank_change>=2), "\n")
cat("max abs rank change:", max(tb$abs_rank_change), "\n\n")
print(as.data.frame(tb), row.names=FALSE, digits=6)
'n countries: 21
n rank changes >=2: 2
max abs rank change: 2
country n raw_mean trimmed_mean raw_rank trimmed_rank
Ethiopia 44 85.4841 85.5178 1 1
United States 10 84.4330 84.9050 2 2
Kenya 25 84.3096 84.5076 3 3
Uganda 36 83.4519 83.4480 4 4
Colombia 183 83.1066 83.2467 5 5
El Salvador 21 83.0529 83.1094 6 6
China 16 82.9275 82.9707 7 8
Costa Rica 51 82.7890 83.0200 8 7
Thailand 32 82.5738 82.6192 9 10
Indonesia 20 82.5655 82.7738 10 9
Peru 10 82.5260 82.4387 11 12
Brazil 132 82.4059 82.5133 12 11
Tanzania, United Republic Of 40 82.3695 82.2719 13 13
Taiwan 75 82.0013 81.9975 14 16
Guatemala 181 81.8466 82.1685 15 14
United States (Hawaii) 73 81.8204 82.0758 16 15
Malawi 11 81.7118 81.7122 17 17
India 14 81.0829 81.3817 18 19
Mexico 236 80.8901 81.2458 19 20
Honduras 52 80.8837 81.5424 20 18
Nicaragua 26 80.4581 80.8941 21 21
abs_rank_change
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) %>% mutate(row_id=as.integer(seq_len(n())-1L))
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
# altitude flags
alt<-dat$altitude_mean_meters; apos<-alt[!is.na(alt)&alt>0]; la<-log10(apos)
q<-quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]; lo_m<-10^(q[[1]]-3*iqr); hi_m<-10^(q[[2]]+3*iqr)
flag_alt<-!is.na(alt)&alt>0&(alt<lo_m|alt>hi_m)
# mv flags
G<-as.matrix(dat[,g7]); mu<-colMeans(G); S<-cov(G); md2<-mahalanobis(G,mu,S)
flag_mv<-md2>qchisq(0.999,df=7)
# influence on grades_defects
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
fr<-dat; fr$d1<-log1p(imp(fr$category_one_defects)); fr$d2<-log1p(imp(fr$category_two_defects)); fr$dq<-log1p(imp(fr$quakers))
f3<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+d1+d2+dq,data=fr)
n<-nobs(f3); p<-length(coef(f3)); hl<-hatvalues(f3)>2*p/n; hc<-cooks.distance(f3)>4/n
infl<-hl&hc
composite<-flag_alt|flag_mv|infl
cat("flag_alt:",sum(flag_alt)," flag_mv:",sum(flag_mv)," infl(both):",sum(infl)," composite:",sum(composite),"\n")
# point estimate
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"])
dp<-100*(cd-cf)/cf
cat(sprintf("coef_full=%.6f coef_drop=%.6f delta_pct=%.6f\n",cf,cd,dp))
cat("n dropped rows:",sum(composite)," n kept:",sum(!composite),"\n")
# bootstrap
set.seed(20260512); B<-500; nn<-nrow(dat); deltas<-numeric(0)
for(b in 1:B){
idx<-sample.int(nn,nn,replace=TRUE)
bd<-dat[idx,]; bc<-composite[idx]
sub<-bd[!bc,]
if(nrow(sub)<10) next
ff<-tryCatch(lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=bd),error=function(e)NULL)
fd<-tryCatch(lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=sub),error=function(e)NULL)
if(is.null(ff)||is.null(fd)) next
a<-unname(coef(ff)["flavor"]); d<-unname(coef(fd)["flavor"])
if(is.na(a)||is.na(d)||!is.finite(a)||!is.finite(d)||a==0) next
deltas<-c(deltas,100*(d-a)/a)
}
ci<-quantile(deltas,c(.025,.975))
cat(sprintf("n_boot_success=%d CI=[%.6f, %.6f]\n",length(deltas),ci[[1]],ci[[2]]))
'flag_alt: 51 flag_mv: 40 infl(both): 35 composite: 105 coef_full=2.085398 coef_drop=1.652728 delta_pct=-20.747640 n dropped rows: 105 n kept: 1233 n_boot_success=500 CI=[-44.328106, 6.247039]
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(stringr)})
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")
allzero <- rowSums(raw[,g10]==0)==10
cat("all-ten-zero rows:", sum(allzero), " | tcp==0 rows:", sum(raw$total_cup_points==0),"\n")
dat <- raw[!allzero,]
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
alt<-dat$altitude_mean_meters; apos<-alt[!is.na(alt)&alt>0]; la<-log10(apos)
q<-quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]; lo_m<-10^(q[[1]]-3*iqr); hi_m<-10^(q[[2]]+3*iqr)
flag_alt<-!is.na(alt)&alt>0&(alt<lo_m|alt>hi_m)
cat(sprintf("fence meters: [%.10f, %.10f]\n", lo_m, hi_m))
G<-as.matrix(dat[,g7]); md2<-mahalanobis(G,colMeans(G),cov(G)); flag_mv<-md2>qchisq(0.999,df=7)
imp<-function(x){x[is.na(x)]<-mean(x,na.rm=TRUE);x}
fr<-dat; fr$d1<-log1p(imp(fr$category_one_defects)); fr$d2<-log1p(imp(fr$category_two_defects)); fr$dq<-log1p(imp(fr$quakers))
f3<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+d1+d2+dq,data=fr)
n<-nobs(f3); p<-length(coef(f3)); hl<-hatvalues(f3)>2*p/n; hc<-cooks.distance(f3)>4/n; infl<-hl&hc
df<-function(x,th){!is.na(x)&x>th}
D<-df(dat$category_one_defects,quantile(dat$category_one_defects,.99,type=7))|df(dat$category_two_defects,quantile(dat$category_two_defects,.99,na.rm=T,type=7))|df(dat$quakers,quantile(dat$quakers,.99,na.rm=T,type=7))
A<-flag_alt; M<-flag_mv; I<-infl
tot<-as.integer(A)+as.integer(M)+as.integer(D)+as.integer(I)
cat("counts by #flags 0..4:", sapply(0:4,function(k)sum(tot==k)),"\n")
cat("alt_only:",sum(A&!M&!I&!D)," mv_only:",sum(M&!A&!I&!D)," infl_only:",sum(I&!A&!M&!D)," defect_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("defect union total:",sum(D),"\n")
'all-ten-zero rows: 1 | tcp==0 rows: 1 fence meters: [357.4462890625, 4923.8166791886] counts by #flags 0..4: 1208 104 21 5 0 alt_only: 43 mv_only: 24 infl_only: 12 defect_only: 25 alt&mv: 4 alt&infl: 5 mv&infl: 13 alt&mv&infl: 1 defect union total: 35
command: create
path: /app/analysis.R
file_text: #!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-matched rewrite
# =====================================================
# Each column is audited with a method that matches its distribution rather
# than a single off-the-shelf summary. See audit_memo.md for the rationale
# behind every block. This script alone reproduces all artifacts in
# /app/outputs from a clean directory.
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)
GRADES7 <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "cupper_points")
GRADES10 <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"uniformity", "clean_cup", "sweetness", "cupper_points")
DEFECT_COLS <- c("category_one_defects", "category_two_defects", "quakers")
# ---------------------------------------------------------------------------
# 1. Load + sentinel drop (a single withdrawn lot has every grade == 0)
# ---------------------------------------------------------------------------
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
all_zero <- rowSums(raw[, GRADES10] == 0, na.rm = TRUE) == length(GRADES10)
dat <- raw[!all_zero, , drop = FALSE]
n_after <- nrow(dat)
# 0-indexed post-drop position, preserving input order
dat$row_id <- as.integer(seq_len(n_after) - 1L)
# ---------------------------------------------------------------------------
# 2. Altitude -- Tukey fence on the log10 scale (k = 3), back-transformed
# ---------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
pos <- !is.na(alt) & alt > 0
log_alt_pos <- log10(alt[pos])
qa <- quantile(log_alt_pos, c(0.25, 0.75), names = FALSE)
iqr_log <- qa[2] - qa[1]
log_lo <- qa[1] - 3 * iqr_log
log_hi <- qa[2] + 3 * iqr_log
alt_lower_m <- 10 ^ log_lo
alt_upper_m <- 10 ^ log_hi
altitude_outlier_flag <- pos & (alt < alt_lower_m | alt > alt_upper_m)
# Recover decimal-displacement unit slips from the raw `altitude` string:
# first numeric token, tested as /10, then /100, then as-is; keep the first
# candidate that lands inside the fence.
first_numeric_token <- function(s) {
as.numeric(str_extract(s, "[0-9]+(?:\\.[0-9]+)?"))
}
altitude_corrected_m <- rep(NA_real_, n_after)
for (i in which(altitude_outlier_flag)) {
tok <- first_numeric_token(dat$altitude[i])
if (is.na(tok)) next
candidates <- c(tok / 10, tok / 100, tok)
inb <- candidates[candidates >= alt_lower_m & candidates <= alt_upper_m]
if (length(inb) > 0) altitude_corrected_m[i] <- inb[1]
}
n_alt_flagged <- sum(altitude_outlier_flag)
n_unit_corrected <- sum(!is.na(altitude_corrected_m))
# elevation used downstream: corrected metres where recovered, else raw mean
elev_used <- ifelse(!is.na(altitude_corrected_m), altitude_corrected_m, alt)
# ---------------------------------------------------------------------------
# 3. Defects -- per-column 99th-percentile upper threshold (mass-at-zero)
# ---------------------------------------------------------------------------
defect_flag_col <- function(x) {
thr <- as.numeric(quantile(x, 0.99, na.rm = TRUE, type = 7))
list(upper_threshold = thr, flag = !is.na(x) & x > thr)
}
d_cat1 <- defect_flag_col(dat$category_one_defects)
d_cat2 <- defect_flag_col(dat$category_two_defects)
d_quak <- defect_flag_col(dat$quakers)
defect_outlier_flag <- d_cat1$flag | d_cat2$flag | d_quak$flag
# ---------------------------------------------------------------------------
# 4. Multivariate grades -- classical Mahalanobis on the 7 primary attributes
# ---------------------------------------------------------------------------
G <- as.matrix(dat[, GRADES7])
cc <- complete.cases(G)
mu <- colMeans(G[cc, , drop = FALSE])
S <- cov(G[cc, , drop = FALSE])
m
…[truncated]Existing file was successfully overwritten at /app/analysis.R
cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1
Error: unexpected 'else' in " else" Execution halted
command: str_replace
path: /app/analysis.R
old_str: selected_model <- if (length(finite_aic) == 0) "grades_only"
else names(finite_aic)[which.min(finite_aic)]
new_str: selected_model <- if (length(finite_aic) == 0) {
"grades_only"
} else {
names(finite_aic)[which.min(finite_aic)]
}The file /app/analysis.R has been edited successfully.
command: str_replace
path: /app/analysis.R
old_str: boot_ci <- if (n_boot > 0) quantile(boot_delta, c(0.025, 0.975), names = FALSE)
else c(NA_real_, NA_real_)
new_str: boot_ci <- if (n_boot > 0) {
quantile(boot_delta, c(0.025, 0.975), names = FALSE)
} else {
c(NA_real_, NA_real_)
}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/
Audit complete: n_input=1339 n_after=1338 selected=grades_defects composite=105 boot=500 ---FILES--- total 225 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:23 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:23 .. -rw-r--r-- 1 h2tagent h2tagent 137149 Jun 20 18:23 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78297 Jun 20 18:23 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6652 Jun 20 18:23 outlier_report.json
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), back-transformed to metres; flagged rows unit-corrected via first-numeric-token /10,/100,as-is",
"defects": "Per-column 99th-percentile upper threshold on zero-inflated counts; row flagged if any of category_one/two/quakers exceeds its threshold",
"multivariate_grades": "Classical Mahalanobis distance (sample mean/covariance) on the 7 primary SCA attributes vs chi-square(0.999, df=7)",
"influence": "Lowest-finite-AIC OLS of 3 specs; on selected fit flag leverage>2p/n AND Cook's D>4/n (both required)",
"country_ranking": "Countries with n>=10: raw mean vs 10% trimmed mean of total_cup_points, compared by rank"
},
"altitude": {
"lower_bound_m": 357.446289062499,
"upper_bound_m": 4923.8166791886,
"n_flagged": 51,
"n_unit_corrected": 7
},
"defects": {
"category_one_defects": {
"upper_threshold": 7.62999999999988,
"n_flagged": 14
},
"category_two_defects": {
"upper_threshold": 26.6299999999999,
"n_flagged": 14
},
"quakers": {
"upper_threshold": 5,
"n_flagged": 10
}
},
"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
…[truncated]cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
js <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
cat("=== CSV structure ===\n")
cat("nrow:", nrow(csv), " (expect 1338)\n")
cat("cols:", paste(names(csv),collapse=", "),"\n")
cat("row_id range:", min(csv$row_id), max(csv$row_id), " monotone:", all(diff(csv$row_id)==1),"\n")
cat("altitude_corrected_m non-NA:", sum(!is.na(csv$altitude_corrected_m))," (expect 7)\n")
cat("corrected values:", paste(sort(csv$altitude_corrected_m[!is.na(csv$altitude_corrected_m)]),collapse=", "),"\n\n")
cat("=== Reload rds & recompute diagnostics ===\n")
fit <- readRDS("outputs/coffee_ols_selected.rds")
cat("class:", class(fit), " nobs:", nobs(fit), " p:", length(coef(fit)),"\n")
cat("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
cat("recomputed n_high_lev:",sum(hl)," n_high_cook:",sum(hc)," both:",sum(hl&hc),"\n")
cat("CSV high_leverage:",sum(csv$high_leverage_flag)," high_cook:",sum(csv$high_cook_flag),"\n")
cat("match leverage:", all(hl==csv$high_leverage_flag)," match cook:", all(hc==csv$high_cook_flag),"\n")
# top10 cook
ord<-order(ck,decreasing=TRUE)[1:10]
cat("recomputed top10 row_ids (0-idx pos):", paste(ord-1,collapse=","),"\n")
cat("JSON top10:", paste(unlist(js$influence$top10_cook_row_ids),collapse=","),"\n")
# composite reconcile
comp <- csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag | (csv$high_leverage_flag & csv$high_cook_flag)
cat("composite match:", all(comp==csv$composite_drop_flag), " sum:",sum(csv$composite_drop_flag),"\n")
'=== CSV structure === nrow: 1338 (expect 1338) 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 row_id range: 0 1337 monotone: TRUE altitude_corrected_m non-NA: 7 (expect 7) corrected values: 518, 518, 518, 1100, 1100, 1901.64, 1901.64 === Reload rds & recompute diagnostics === class: lm nobs: 1338 p: 11 terms: (Intercept), aroma, flavor, aftertaste, acidity, body, balance, cupper_points, d_cat1, d_cat2, d_quak recomputed n_high_lev: 110 n_high_cook: 66 both: 35 CSV high_leverage: 110 high_cook: 66 match leverage: TRUE match cook: TRUE recomputed top10 row_ids (0-idx pos): 1309,1308,1333,1295,1306,1260,1307,1302,1291,1259 JSON top10: 1309,1308,1333,1295,1306,1260,1307,1302,1291,1259 composite match: TRUE sum: 105
cd /app && Rscript -e '
suppressPackageStartupMessages(library(jsonlite))
js <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
errs <- character(0)
chk <- function(cond,msg) if(!isTRUE(cond)) errs<<-c(errs,msg)
need_int <- function(x,nm){ chk(is.numeric(x)&&length(x)==1&&x==round(x), paste(nm,"not int")) }
need_num <- function(x,nm){ chk(is.numeric(x)&&length(x)==1, paste(nm,"not num")) }
need_str <- function(x,nm){ chk(is.character(x)&&length(x)==1, paste(nm,"not str")) }
need_int(js$n_input,"n_input"); need_int(js$n_after_zero_score_drop,"n_after")
for(k in c("altitude","defects","multivariate_grades","influence","country_ranking")) need_str(js$methods[[k]],paste0("methods.",k))
need_num(js$altitude$lower_bound_m,"alt.lower"); need_num(js$altitude$upper_bound_m,"alt.upper")
need_int(js$altitude$n_flagged,"alt.n_flagged"); need_int(js$altitude$n_unit_corrected,"alt.n_unit")
for(c in c("category_one_defects","category_two_defects","quakers")){ need_num(js$defects[[c]]$upper_threshold,paste0("def.",c,".thr")); need_int(js$defects[[c]]$n_flagged,paste0("def.",c,".n")) }
chk(length(js$multivariate_grades$columns_used)==7,"mv cols!=7"); need_num(js$multivariate_grades$threshold,"mv.thr"); need_int(js$multivariate_grades$n_flagged,"mv.n")
need_str(js$influence$selected_model,"infl.sel"); for(k in c("aic_grades_only","aic_grades_altitude","aic_grades_defects")) need_num(js$influence[[k]],paste0("infl.",k))
need_num(js$influence$leverage_threshold,"infl.lev"); need_num(js$influence$cook_threshold,"infl.cook")
need_int(js$influence$n_high_leverage,"infl.nlev"); need_int(js$influence$n_high_cook,"infl.ncook"); need_int(js$influence$n_both,"infl.nboth")
chk(length(js$influence$top10_cook_row_ids)==10,"top10 !=10"); for(v in js$influence$top10_cook_row_ids) need_int(v,"top10 elem")
need_int(js$country_ranking$n_countries_evaluated,"cr.n"); need_int(js$country_ranking$n_rank_changes_ge_2,"cr.ge2"); need_int(js$country_ranking$max_abs_rank_change,"cr.max")
pc<-js$country_ranking$per_country; chk(length(pc)==js$country_ranking$n_countries_evaluated,"per_country len")
for(e in pc){ need_str(e$country,"pc.country"); need_int(e$n,"pc.n"); need_num(e$raw_mean,"pc.raw"); need_num(e$trimmed_mean,"pc.trim"); need_int(e$raw_rank,"pc.rr"); need_int(e$trimmed_rank,"pc.tr") }
for(k in c("coef_flavor_full","coef_flavor_dropped","delta_pct","delta_pct_ci_lower","delta_pct_ci_upper")) need_num(js$sensitivity[[k]],paste0("sens.",k))
need_int(js$sensitivity$n_bootstrap_replicates,"sens.nboot")
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_int(js$flag_interactions[[k]],paste0("fi.",k))
if(length(errs)==0) cat("ALL SCHEMA CHECKS PASSED\n") else { cat("ERRORS:\n"); print(errs) }
'ALL SCHEMA CHECKS PASSED
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw[!(rowSums(raw[,c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")]==0)==10),]
js <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
# defect consistency
for(c in c("category_one_defects","category_two_defects","quakers")){
thr<-js$defects[[c]]$upper_threshold; n<-js$defects[[c]]$n_flagged
act<-sum(!is.na(dat[[c]]) & dat[[c]]>thr)
cat(sprintf("%-22s thr=%.4f json_n=%d recompute=%d match=%s\n",c,thr,n,act,n==act))
}
csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
cat("\naltitude flag CSV vs JSON:",sum(csv$altitude_outlier_flag),"==",js$altitude$n_flagged,"\n")
cat("mv flag CSV vs JSON:",sum(csv$multivariate_grade_outlier_flag),"==",js$multivariate_grades$n_flagged,"\n")
cat("defect flag CSV total:",sum(csv$defect_outlier_flag),"\n")
fi<-js$flag_interactions
cat("\nflag-count partition sums to n_after:",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,"(expect 1338)\n")
# one-flag breakdown should be <= n_with_one_flag and equal sum of *_only
cat("sum of *_only:",fi$n_alt_only+fi$n_mv_only+fi$n_influence_only+fi$n_defect_only,"== n_with_one_flag:",fi$n_with_one_flag,"\n")
'
echo "=== Determinism: re-run and diff JSON ==="
cd /app && cp outputs/outlier_report.json /tmp/run1.json && rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1 && diff -q /tmp/run1.json outputs/outlier_report.json && echo "JSON identical across runs"category_one_defects thr=7.6300 json_n=14 recompute=14 match=TRUE category_two_defects thr=26.6300 json_n=14 recompute=14 match=TRUE quakers thr=5.0000 json_n=10 recompute=10 match=TRUE altitude flag CSV vs JSON: 51 == 51 mv flag CSV vs JSON: 40 == 40 defect flag CSV total: 35 flag-count partition sums to n_after: 1338 (expect 1338) sum of *_only: 104 == n_with_one_flag: 104 === Determinism: re-run and diff JSON === JSON identical across runs
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw[raw$total_cup_points>0,]
for(c in c("uniformity","clean_cup","sweetness")){
x<-dat[[c]]; cat(sprintf("%-11s pct==10: %.1f%% min=%.1f\n",c,100*mean(x==10,na.rm=TRUE),min(x,na.rm=TRUE)))
}
'
echo "=== CSV + rds determinism ==="
cd /app && cp outputs/outlier_flags.csv /tmp/f1.csv && md5sum outputs/coffee_ols_selected.rds > /tmp/r1.md5
rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1
diff -q /tmp/f1.csv outputs/outlier_flags.csv && echo "CSV identical"
Rscript -e 'f<-readRDS("outputs/coffee_ols_selected.rds"); cat("reload OK, coef flavor:",round(coef(f)["flavor"],5),"\n")'uniformity pct==10: 86.1% min=6.0 clean_cup pct==10: 91.1% min=0.0 sweetness pct==10: 91.0% min=1.3 === CSV + rds determinism === CSV identical reload OK, coef flavor: 2.08014
command: create path: /app/outputs/audit_memo.md file_text: # Coffee Quality Outlier Audit , Method Memo The draft applied one off-the-shelf summary to every column. Each block below says why that default is wrong for the column's actual shape and what replaced it. All figures come from re-running `/app/analysis.R`. ## Data prep The draft filtered `total_cup_points > 0`. That happens to remove the single withdrawn lot here, but it conflates "total score is zero" with the real signal , *every* grade recorded as zero. I drop the one lot whose ten grade columns are all zero (a withdrawn submission), which is the documented sentinel and is robust to any future low-but-valid lot. Input = 1339, post-drop = 1338. `row_id` is reassigned 0..1337 in post-drop input order so every downstream artifact shares one index. ## Altitude The draft put a Tukey *k*=3 fence on raw metres. `altitude_mean_meters` is strongly right-skewed (median 1310 m, max 190164 m), so a symmetric IQR fence on the raw scale is dominated by the large tail and mis-describes what is a multiplicative quantity; it also offers no path to fixing unit slips. I take `log10` of positive altitudes, build the *k*=3 IQR fence there, and back-transform to **[357 m, 4924 m]**; 51 rows fall outside. I then try to recover decimal-displacement errors from the raw `altitude` text: the first numeric token tested as ÷10, ÷100, then as-is, keeping the first candidate inside the fence. That recovers **7** lots (e.g. `11000 metros`→1100, `190164`→1901.64, `1100.00 mosl`→1100, `518` ft→518); the rest stay `NA`. The draft never attempted recovery (`n_unit_corrected = 0`). ## Defects The draft used raw Tukey *k*=3 fences. Because 85% of `category_one_defects` and 93% of `quakers` are exactly 0, Q1=Q3=0, the IQR collapses to 0, and the fence degenerates to "> 0" , flagging 202 and 94 rows (15% / 7%). That is noise, not extremeness. These are zero-inflated counts with a sparse upper tail, so I set a per-column **99th-percentile** upper threshold and flag counts strictly above it: thresholds 7.63 / 26.63 / 5 giving 14 / 14 / 10 rows (~1%). A lot is a defect outlier if any column trips (35 rows). This targets genuinely extreme counts while ignoring the mass at zero. ## Multivariate grades The draft ran Mahalanobis on all ten grade columns. Three of them , `uniformity`, `clean_cup`, `sweetness` , are near-constant (86–91% sit exactly at 10.0 with a few discrete drops). Near-constant columns make the covariance ill-conditioned (the draft needed `tol = 1e-30`), and their rare low values dominate the distance, so a chi-square(df=10) reference is invalid. I restrict to the seven primary SCA attributes (well-conditioned covariance), use a **classical** mean/covariance Mahalanobis distance (not MCD, per spec), and threshold at chi-square(0.999, df=7) = 24.32. That flags **40** joint outliers. ## Influence The draft used `abs(rstandard) > 2` and reported it simultaneously as leverage, Cook's D, and "both" , three distinct diagnostics reduced to one large-residual rule that flags ~5% by construction and never measures leverage or actual influence. It also considered only one model. I fit three OLS specs for `total_cup_points` and pick the lowest finite AIC: grades_only (4459.4), grades_altitude (4461.3), grades_defects (**4418.7**, selected) , the `log1p` defect terms carry real signal. On the selected fit I compute hat-value leverage (> 2p/n = 0.0164 → 110 rows) and Cook's D (> 4/n = 0.00299 → 66 rows) **separately**, and define influence as the intersection (**35** rows). The ten largest Cook's D positions are reported for follow-up. ## Country ranking The draft reported raw means and hard-coded the trimmed mean equal to the raw mean (every rank change = 0) , no robustness check at all. For the 21 countries with ≥10 lots I compute a 10% trimmed mean beside the raw mean and rank by each. Two countries move ≥2 positions (max shift 2): Taiwan falls 14→16 and Honduras rises 20→18 once low-scoring lo …[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 '^## ' outputs/audit_memo.md && echo "--- Expected 8 ---" && printf '%s\n' "## Data prep" "## Altitude" "## Defects" "## Multivariate grades" "## Influence" "## Country ranking" "## Sensitivity" "## Composite flags" | diff - <(grep '^## ' outputs/audit_memo.md) && echo "HEADINGS MATCH EXACTLY"
Word count: 907 --- Headings --- ## Data prep ## Altitude ## Defects ## Multivariate grades ## Influence ## Country ranking ## Sensitivity ## Composite flags --- Expected 8 --- HEADINGS MATCH EXACTLY
command: str_replace path: /app/outputs/audit_memo.md old_str: The draft applied one off-the-shelf summary to every column. Each block below says why that default is wrong for the column's actual shape and what replaced it. All figures come from re-running `/app/analysis.R`. new_str: The draft applied one off-the-shelf summary to every column. Each block says why that default is wrong for the column's shape and what replaced it.
The file /app/outputs/audit_memo.md has been edited successfully.
command: str_replace path: /app/outputs/audit_memo.md old_str: signal , *every* grade recorded as zero. I drop the one lot whose ten grade columns are all zero (a withdrawn submission), which is the documented sentinel and is robust to any future low-but-valid lot. Input = 1339, new_str: signal , *every* grade recorded as zero. I drop the one lot whose ten grade columns are all zero (the withdrawn sentinel), robust to any future low-but-valid lot. Input = 1339,
The file /app/outputs/audit_memo.md has been edited successfully.
cd /app && echo "Word count: $(wc -w < outputs/audit_memo.md)"
Word count: 892
cd /app && grep -n ')---"' outputs/audit_memo.md || echo "delimiter )---\" not present, safe to use r\"---( )---\"" echo "--- check for )----- sequence ---" && grep -n ')-----' outputs/audit_memo.md || echo "safe"
delimiter )---" not present, safe to use r"---( )---" --- check for )----- sequence --- safe
cd /app && cat outputs/audit_memo.md
# Coffee Quality Outlier Audit , Method Memo The draft applied one off-the-shelf summary to every column. Each block says why that default is wrong for the column's shape and what replaced it. ## Data prep The draft filtered `total_cup_points > 0`. That happens to remove the single withdrawn lot here, but it conflates "total score is zero" with the real signal , *every* grade recorded as zero. I drop the one lot whose ten grade columns are all zero (the withdrawn sentinel), robust to any future low-but-valid lot. Input = 1339, post-drop = 1338. `row_id` is reassigned 0..1337 in post-drop input order so every downstream artifact shares one index. ## Altitude The draft put a Tukey *k*=3 fence on raw metres. `altitude_mean_meters` is strongly right-skewed (median 1310 m, max 190164 m), so a symmetric IQR fence on the raw scale is dominated by the large tail and mis-describes what is a multiplicative quantity; it also offers no path to fixing unit slips. I take `log10` of positive altitudes, build the *k*=3 IQR fence there, and back-transform to **[357 m, 4924 m]**; 51 rows fall outside. I then try to recover decimal-displacement errors from the raw `altitude` text: the first numeric token tested as ÷10, ÷100, then as-is, keeping the first candidate inside the fence. That recovers **7** lots (e.g. `11000 metros`→1100, `190164`→1901.64, `1100.00 mosl`→1100, `518` ft→518); the rest stay `NA`. The draft never attempted recovery (`n_unit_corrected = 0`). ## Defects The draft used raw Tukey *k*=3 fences. Because 85% of `category_one_defects` and 93% of `quakers` are exactly 0, Q1=Q3=0, the IQR collapses to 0, and the fence degenerates to "> 0" , flagging 202 and 94 rows (15% / 7%). That is noise, not extremeness. These are zero-inflated counts with a sparse upper tail, so I set a per-column **99th-percentile** upper threshold and flag counts strictly above it: thresholds 7.63 / 26.63 / 5 giving 14 / 14 / 10 rows (~1%). A lot is a defect outlier if any column trips (35 rows). This targets genuinely extreme counts while ignoring the mass at zero. ## Multivariate grades The draft ran Mahalanobis on all ten grade columns. Three of them , `uniformity`, `clean_cup`, `sweetness` , are near-constant (86–91% sit exactly at 10.0 with a few discrete drops). Near-constant columns make the covariance ill-conditioned (the draft needed `tol = 1e-30`), and their rare low values dominate the distance, so a chi-square(df=10) reference is invalid. I restrict to the seven primary SCA attributes (well-conditioned covariance), use a **classical** mean/covariance Mahalanobis distance (not MCD, per spec), and threshold at chi-square(0.999, df=7) = 24.32. That flags **40** joint outliers. ## Influence The draft used `abs(rstandard) > 2` and reported it simultaneously as leverage, Cook's D, and "both" , three distinct diagnostics reduced to one large-residual rule that flags ~5% by construction and never measures leverage or actual influence. It also considered only one model. I fit three OLS specs for `total_cup_points` and pick the lowest finite AIC: grades_only (4459.4), grades_altitude (4461.3), grades_defects (**4418.7**, selected) , the `log1p` defect terms carry real signal. On the selected fit I compute hat-value leverage (> 2p/n = 0.0164 → 110 rows) and Cook's D (> 4/n = 0.00299 → 66 rows) **separately**, and define influence as the intersection (**35** rows). The ten largest Cook's D positions are reported for follow-up. ## Country ranking The draft reported raw means and hard-coded the trimmed mean equal to the raw mean (every rank change = 0) , no robustness check at all. For the 21 countries with ≥10 lots I compute a 10% trimmed mean beside the raw mean and rank by each. Two countries move ≥2 positions (max shift 2): Taiwan falls 14→16 and Honduras rises 20→18 once low-scoring lots are trimmed. The ranking is mostly stable, but these shifts show a few national means are pulled by low outliers , exactly what the draft could not detect. ## Se …[truncated]
cd /app && Rscript -e 'cat("locale:", Sys.getlocale("LC_CTYPE"), "\n"); cat("native enc:", l10n_info()$`UTF-8`, "\n")'locale: en_US.UTF-8 native enc: TRUE
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)
cat(sprintf("Audit complete: n_input=%d n_after=%d selected=%s composite=%d boot=%d\n",
n_input, n_after, selected_model, sum(composite_drop_flag), n_boot))
new_str: write_json(results, file.path(OUTPUT_DIR, "outlier_report.json"),
auto_unbox = FALSE, pretty = TRUE, na = "null", digits = NA)
# ---------------------------------------------------------------------------
# Method memo (<=900 words). Written by this script so the artifact set is
# fully reproducible from analysis.R alone.
# ---------------------------------------------------------------------------
memo_text <- enc2utf8(r"---(# Coffee Quality Outlier Audit , Method Memo
The draft applied one off-the-shelf summary to every column. Each block says
why that default is wrong for the column's shape and what replaced it.
## Data prep
The draft filtered `total_cup_points > 0`. That happens to remove the single
withdrawn lot here, but it conflates "total score is zero" with the real
signal , *every* grade recorded as zero. I drop the one lot whose ten grade
columns are all zero (the withdrawn sentinel), robust to any future
low-but-valid lot. Input = 1339,
post-drop = 1338. `row_id` is reassigned 0..1337 in post-drop input order so
every downstream artifact shares one index.
## Altitude
The draft put a Tukey *k*=3 fence on raw metres. `altitude_mean_meters` is
strongly right-skewed (median 1310 m, max 190164 m), so a symmetric IQR fence
on the raw scale is dominated by the large tail and mis-describes what is a
multiplicative quantity; it also offers no path to fixing unit slips. I take
`log10` of positive altitudes, build the *k*=3 IQR fence there, and
back-transform to **[357 m, 4924 m]**; 51 rows fall outside. I then try to
recover decimal-displacement errors from the raw `altitude` text: the first
numeric token tested as ÷10, ÷100, then as-is, keeping the first candidate
inside the fence. That recovers **7** lots (e.g. `11000 metros`→1100,
`190164`→1901.64, `1100.00 mosl`→1100, `518` ft→518); the rest stay `NA`. The
draft never attempted recovery (`n_unit_corrected = 0`).
## Defects
The draft used raw Tukey *k*=3 fences. Because 85% of `category_one_defects`
and 93% of `quakers` are exactly 0, Q1=Q3=0, the IQR collapses to 0, and the
fence degenerates to "> 0" , flagging 202 and 94 rows (15% / 7%). That is
noise, not extremeness. These are zero-inflated counts with a sparse upper
tail, so I set a per-column **99th-percentile** upper threshold and flag counts
strictly above it: thresholds 7.63 / 26.63 / 5 giving 14 / 14 / 10 rows (~1%).
A lot is a defect outlier if any column trips (35 rows). This targets
genuinely extreme counts while ignoring the mass at zero.
## Multivariate grades
The draft ran Mahalanobis on all ten grade columns. Three of them ,
`uniformity`, `clean_cup`, `sweetness` , are near-constant (86–91% sit exactly
at 10.0 with a few discrete drops). Near-constant columns make the covariance
ill-conditioned (the draft needed `tol = 1e-30`), and their rare low values
dominate the distance, so a chi-square(df=10) reference is invalid. I restrict
to the seven primary SCA attributes (well-conditioned covariance), use a
**classical** mean/covariance Mahalanobis distance (not MCD, per spec), and
threshold at chi-square(0.999, df=7) = 24.32. That flags **40** joint outliers.
## Influence
The draft used `abs(rstandard) > 2` and reported it simultaneously as leverage,
Cook's D, and "both" , three distinct diagnostics reduced to one large-residual
rule that flags ~5% by construction and never measures leverage or actual
influence. It also considered only one model. I fit three OLS specs for
`total_cup_points` and pick the lowest finite AIC: grades_only (4459.4),
grades_altitude (4461.3), grades_defects (**4418.7**, selected) , the `log1p`
defect te
…[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 word count---" && wc -w < outputs/audit_memo.md && echo "---MEMO headings---" && grep -c '^## ' outputs/audit_memo.md
Audit complete: n_input=1339 n_after=1338 selected=grades_defects composite=105 boot=500 ---FILES--- total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:26 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:26 .. -rw-r--r-- 1 h2tagent h2tagent 5635 Jun 20 18:26 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137149 Jun 20 18:26 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78297 Jun 20 18:26 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6652 Jun 20 18:26 outlier_report.json ---MEMO word count--- 892 ---MEMO headings--- 8
cd /app && echo "--- heading exact match ---" && printf '%s\n' "## Data prep" "## Altitude" "## Defects" "## Multivariate grades" "## Influence" "## Country ranking" "## Sensitivity" "## Composite flags" | diff - <(grep '^## ' outputs/audit_memo.md) && echo "HEADINGS MATCH" && echo "--- UTF-8 sample lines ---" && grep -n '÷10\|→1100\|≥10\|∧ influence\|−20.7' outputs/audit_memo.md && echo "--- file encoding ---" && file outputs/audit_memo.md
--- heading exact match --- HEADINGS MATCH --- UTF-8 sample lines --- 23:numeric token tested as ÷10, ÷100, then as-is, keeping the first candidate 24:inside the fence. That recovers **7** lots (e.g. `11000 metros`→1100, 25:`190164`→1901.64, `1100.00 mosl`→1100, `518` ft→518); the rest stay `NA`. The 63:with ≥10 lots I compute a 10% trimmed mean beside the raw mean and rank by 74:`flavor` coefficient 2.085 → 1.653, a **−20.7%** change. I bootstrap that 87:is multivariate ∧ influence (13) , lots with extreme grade profiles also move 88:the fit , while altitude ∧ multivariate (4) and altitude ∧ influence (5) are --- file encoding --- outputs/audit_memo.md: Unicode text, UTF-8 text
cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1
cp outputs/outlier_report.json /tmp/a.json; cp outputs/outlier_flags.csv /tmp/a.csv; cp outputs/audit_memo.md /tmp/a.md
rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1
echo "=== Determinism across runs ==="
diff -q /tmp/a.json outputs/outlier_report.json && echo "JSON identical"
diff -q /tmp/a.csv outputs/outlier_flags.csv && echo "CSV identical"
diff -q /tmp/a.md outputs/audit_memo.md && echo "MEMO identical"
echo ""
echo "=== Full grader-style reconciliation (fresh session) ==="
Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite)})
fit <- readRDS("outputs/coffee_ols_selected.rds")
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
js <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
n<-nobs(fit); p<-length(coef(fit))
stopifnot(js$influence$selected_model=="grades_defects")
stopifnot(abs(js$influence$leverage_threshold-2*p/n)<1e-12)
stopifnot(abs(js$influence$cook_threshold-4/n)<1e-12)
hl<-hatvalues(fit)>2*p/n; hc<-cooks.distance(fit)>4/n
stopifnot(all(hl==csv$high_leverage_flag), all(hc==csv$high_cook_flag))
stopifnot(sum(hl)==js$influence$n_high_leverage, sum(hc)==js$influence$n_high_cook, sum(hl&hc)==js$influence$n_both)
ord<-as.integer(order(cooks.distance(fit),decreasing=TRUE)[1:10]-1)
stopifnot(identical(ord, as.integer(unlist(js$influence$top10_cook_row_ids))))
comp <- csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag | (csv$high_leverage_flag & csv$high_cook_flag)
stopifnot(all(comp==csv$composite_drop_flag))
stopifnot(sum(!is.na(csv$altitude_corrected_m))==js$altitude$n_unit_corrected)
# CSV <-> input order check
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")
din<-raw[!(rowSums(raw[,g10]==0,na.rm=TRUE)==10),]
stopifnot(nrow(csv)==nrow(din), all(csv$total_cup_points==din$total_cup_points))
cat("ALL RECONCILIATION CHECKS PASSED (n=",nrow(csv),", p=",p,")\n",sep="")
'=== Determinism across runs === JSON identical CSV identical MEMO identical === Full grader-style reconciliation (fresh session) === ALL RECONCILIATION CHECKS PASSED (n=1338, p=11)
cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(jsonlite)})
js <- fromJSON("outputs/outlier_report.json", simplifyVector=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),]
# altitude bound round-trip
alt<-dat$altitude_mean_meters; la<-log10(alt[!is.na(alt)&alt>0]); q<-quantile(la,c(.25,.75),names=FALSE); iqr<-q[2]-q[1]
stopifnot(abs(js$altitude$lower_bound_m-10^(q[1]-3*iqr))<1e-6, abs(js$altitude$upper_bound_m-10^(q[2]+3*iqr))<1e-6)
cat("altitude bounds round-trip OK\n")
# independent country recompute
tb <- dat %>% filter(!is.na(country_of_origin)) %>% group_by(country=country_of_origin) %>%
summarise(n=n(), rm=mean(total_cup_points), tm=mean(total_cup_points,trim=0.1), .groups="drop") %>%
filter(n>=10) %>% mutate(rr=rank(-rm,ties.method="min"), tr=rank(-tm,ties.method="min"))
stopifnot(nrow(tb)==js$country_ranking$n_countries_evaluated)
stopifnot(sum(abs(tb$rr-tb$tr)>=2)==js$country_ranking$n_rank_changes_ge_2)
stopifnot(max(abs(tb$rr-tb$tr))==js$country_ranking$max_abs_rank_change)
# check each per_country entry against recompute
pc<-js$country_ranking$per_country
for(e in pc){ row<-tb[tb$country==e$country,]; stopifnot(nrow(row)==1,
row$n==e$n, abs(row$rm-e$raw_mean)<1e-9, abs(row$tm-e$trimmed_mean)<1e-9,
row$rr==e$raw_rank, row$tr==e$trimmed_rank) }
cat("country ranking fully reconciled (", nrow(tb), "countries )\n")
# sensitivity point estimate sign/containment
s<-js$sensitivity
stopifnot(s$delta_pct_ci_lower < s$delta_pct_ci_upper)
stopifnot(s$delta_pct >= s$delta_pct_ci_lower, s$delta_pct <= s$delta_pct_ci_upper)
stopifnot(s$n_bootstrap_replicates==500)
cat(sprintf("sensitivity OK: delta=%.3f in CI [%.3f, %.3f], reps=%d\n",
s$delta_pct, s$delta_pct_ci_lower, s$delta_pct_ci_upper, s$n_bootstrap_replicates))
cat("\nALL INDEPENDENT CROSS-CHECKS PASSED\n")
'altitude bounds round-trip OK country ranking fully reconciled ( 21 countries ) sensitivity OK: delta=-20.748 in CI [-44.328, 6.247], reps=500 ALL INDEPENDENT CROSS-CHECKS PASSED
cd /app && rm -rf outputs && mkdir outputs && Rscript -e 'options(warn=1); source("analysis.R")' 2>&1 | grep -iE "warning|error" || echo "No warnings or errors emitted."No warnings or errors emitted.
cd /app && sed -n '1,130p' analysis.R
#!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-matched rewrite
# =====================================================
# Each column is audited with a method that matches its distribution rather
# than a single off-the-shelf summary. See audit_memo.md for the rationale
# behind every block. This script alone reproduces all artifacts in
# /app/outputs from a clean directory.
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)
GRADES7 <- c("aroma", "flavor", "aftertaste", "acidity",
"body", "balance", "cupper_points")
GRADES10 <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance",
"uniformity", "clean_cup", "sweetness", "cupper_points")
DEFECT_COLS <- c("category_one_defects", "category_two_defects", "quakers")
# ---------------------------------------------------------------------------
# 1. Load + sentinel drop (a single withdrawn lot has every grade == 0)
# ---------------------------------------------------------------------------
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)
all_zero <- rowSums(raw[, GRADES10] == 0, na.rm = TRUE) == length(GRADES10)
dat <- raw[!all_zero, , drop = FALSE]
n_after <- nrow(dat)
# 0-indexed post-drop position, preserving input order
dat$row_id <- as.integer(seq_len(n_after) - 1L)
# ---------------------------------------------------------------------------
# 2. Altitude -- Tukey fence on the log10 scale (k = 3), back-transformed
# ---------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
pos <- !is.na(alt) & alt > 0
log_alt_pos <- log10(alt[pos])
qa <- quantile(log_alt_pos, c(0.25, 0.75), names = FALSE)
iqr_log <- qa[2] - qa[1]
log_lo <- qa[1] - 3 * iqr_log
log_hi <- qa[2] + 3 * iqr_log
alt_lower_m <- 10 ^ log_lo
alt_upper_m <- 10 ^ log_hi
altitude_outlier_flag <- pos & (alt < alt_lower_m | alt > alt_upper_m)
# Recover decimal-displacement unit slips from the raw `altitude` string:
# first numeric token, tested as /10, then /100, then as-is; keep the first
# candidate that lands inside the fence.
first_numeric_token <- function(s) {
as.numeric(str_extract(s, "[0-9]+(?:\\.[0-9]+)?"))
}
altitude_corrected_m <- rep(NA_real_, n_after)
for (i in which(altitude_outlier_flag)) {
tok <- first_numeric_token(dat$altitude[i])
if (is.na(tok)) next
candidates <- c(tok / 10, tok / 100, tok)
inb <- candidates[candidates >= alt_lower_m & candidates <= alt_upper_m]
if (length(inb) > 0) altitude_corrected_m[i] <- inb[1]
}
n_alt_flagged <- sum(altitude_outlier_flag)
n_unit_corrected <- sum(!is.na(altitude_corrected_m))
# elevation used downstream: corrected metres where recovered, else raw mean
elev_used <- ifelse(!is.na(altitude_corrected_m), altitude_corrected_m, alt)
# ---------------------------------------------------------------------------
# 3. Defects -- per-column 99th-percentile upper threshold (mass-at-zero)
# ---------------------------------------------------------------------------
defect_flag_col <- function(x) {
thr <- as.numeric(quantile(x, 0.99, na.rm = TRUE, type = 7))
list(upper_threshold = thr, flag = !is.na(x) & x > thr)
}
d_cat1 <- defect_flag_col(dat$category_one_defects)
d_cat2 <- defect_flag_col(dat$category_two_defects)
d_quak <- defect_flag_col(dat$quakers)
defect_outlier_flag <- d_cat1$flag | d_cat2$flag | d_quak$flag
# ---------------------------------------------------------------------------
# 4. Multivariate grades -- classical Mahalanobis on the 7 primary attributes
# ---------------------------------------------------------------------------
G <- as.matrix(dat[, GRADES7])
cc <- complete.cases(G)
mu <- colMeans(G[cc, , drop = FALSE])
S <- cov(G[cc, , drop = FALSE])
md2 <- rep(NA_real_, n_after)
md2[cc] <- mahalanob
…[truncated]cd /app && sed -n '130,230p' analysis.R
fit_sel <- switch(selected_model,
grades_only = f_grades,
grades_altitude = f_altitude,
grades_defects = f_defects)
n_fit <- nobs(fit_sel)
p_fit <- length(coef(fit_sel))
lev_thr <- 2 * p_fit / n_fit
cook_thr <- 4 / n_fit
lev <- hatvalues(fit_sel)
cook <- cooks.distance(fit_sel)
high_leverage_local <- lev > lev_thr
high_cook_local <- cook > cook_thr
# map local (complete-case) diagnostics back to all surviving rows
high_leverage_flag <- rep(FALSE, n_after); high_leverage_flag[cc] <- high_leverage_local
high_cook_flag <- rep(FALSE, n_after); high_cook_flag[cc] <- high_cook_local
influence_flag <- high_leverage_flag & high_cook_flag # BOTH criteria
# ten largest Cook's D, reported as 0-indexed post-drop row_ids
top_local <- order(cook, decreasing = TRUE)[seq_len(min(10, length(cook)))]
top10_cook_row_ids <- as.integer(reg$row_id[top_local])
# ---------------------------------------------------------------------------
# 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.1),
.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_rank_change <- if (n_countries > 0) max(country_tbl$abs_rank_change) else 0L
# ---------------------------------------------------------------------------
# 7. Composite drop + flavor-coefficient sensitivity + bootstrap CI
# ---------------------------------------------------------------------------
composite_drop_flag <- altitude_outlier_flag |
multivariate_grade_outlier_flag | influence_flag
grades_formula <- total_cup_points ~ aroma + flavor + aftertaste + acidity +
body + balance + cupper_points
flavor_coef <- function(fit) unname(coef(fit)["flavor"])
fit_full <- lm(grades_formula, data = dat)
fit_dropped <- lm(grades_formula, data = dat[!composite_drop_flag, , drop = FALSE])
coef_full <- flavor_coef(fit_full)
coef_dropped <- flavor_coef(fit_dropped)
delta_pct <- 100 * (coef_dropped - coef_full) / coef_full
set.seed(20260512)
B <- 500
nn <- n_after
boot_delta <- numeric(0)
for (b in seq_len(B)) {
idx <- sample.int(nn, nn, replace = TRUE)
bd <- dat[idx, , drop = FALSE]
bflag <- composite_drop_flag[idx] # carry ORIGINAL flags
sub <- bd[!bflag, , drop = FALSE]
if (nrow(sub) < 10) next
ff <- tryCatch(lm(grades_formula, data = bd), error = function(e) NULL)
fd <- tryCatch(lm(grades_formula, data = sub), error = function(e) NULL)
if (is.null(ff) || is.null(fd)) next
cf <- flavor_coef(ff); cd <- flavor_coef(fd)
if (!is.finite(cf) || !is.finite(cd) || cf == 0) next
boot_delta <- c(boot_delta, 100 * (cd - cf) / cf)
}
n_boot <- length(boot_delta)
boot_ci <- if (n_boot > 0) {
quantile(boot_delta, c(0.025, 0.975), names = FALSE)
} else {
c(NA_real_, NA_real_)
}
# ---------------------------------------------------------------------------
# 8. Flag interactions (influence = leverage AND Cook)
# ---------------------------------------------------------------------------
A <- altitude_outlier_flag
M <- multivariate_grade_outlier_flag
D <- defect_outlier_flag
I <- influence_flag
total_flags <- as.integer(A) + as.integer(M) + as.integer(D) + as.integer(I)
# ---------------------------------------------------------------------------
# Write artifacts
# ---------------------------------------------------------------------------
saveRDS(fit_sel, file.path(OUT
…[truncated]cd /tmp && rm -rf /app/outputs && mkdir -p /app/outputs && Rscript /app/analysis.R 2>&1 && echo "--- final outputs ---" && ls -la /app/outputs/
Audit complete: n_input=1339 n_after=1338 selected=grades_defects composite=105 boot=500 --- final outputs --- total 231 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:28 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:27 .. -rw-r--r-- 1 h2tagent h2tagent 5635 Jun 20 18:28 audit_memo.md -rw-r--r-- 1 h2tagent h2tagent 137149 Jun 20 18:28 coffee_ols_selected.rds -rw-r--r-- 1 h2tagent h2tagent 78297 Jun 20 18:28 outlier_flags.csv -rw-r--r-- 1 h2tagent h2tagent 6652 Jun 20 18:28 outlier_report.json
cd /app && head -5 outputs/outlier_flags.csv && echo "..." && echo "--- rows WITH a unit correction ---" && awk -F, 'NR==1 || $5!="NA"' outputs/outlier_flags.csv | head -9
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 2,Guatemala,89.75,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 3,Ethiopia,89,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE ... --- rows WITH a unit correction --- 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 31,"Tanzania, United Republic Of",86.5,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 109,"Tanzania, United Republic Of",84.92,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 201,"Tanzania, United Republic Of",84.25,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 202,"Tanzania, United Republic Of",84.25,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 221,"Tanzania, United Republic Of",84.17,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 357,"Tanzania, United Republic Of",83.5,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 385,"Tanzania, United Republic Of",83.42,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE 402,"Tanzania, United Republic Of",83.33,FALSE,NA,FALSE,FALSE,FALSE,FALSE,FALSE
cd /app && Rscript -e '
suppressPackageStartupMessages(library(readr))
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
cr <- csv[!is.na(csv$altitude_corrected_m), c("row_id","country_of_origin","altitude_outlier_flag","altitude_corrected_m","composite_drop_flag")]
cat("Rows with a unit correction (", nrow(cr), "):\n", sep="")
print(as.data.frame(cr), row.names=FALSE)
cat("\nAll corrected rows are altitude-flagged:", all(cr$altitude_outlier_flag), "\n")
cat("altitude_corrected_m is NA for all non-corrected rows:", all(is.na(csv$altitude_corrected_m[!csv$row_id %in% cr$row_id])), "\n")
cat("\nColumn classes:\n"); print(sapply(csv, class))
'Rows with a unit correction (7):
row_id country_of_origin altitude_outlier_flag altitude_corrected_m
543 Brazil TRUE 1100.00
584 Guatemala TRUE 518.00
677 Guatemala TRUE 518.00
833 Guatemala TRUE 518.00
896 Guatemala TRUE 1901.64
1040 Nicaragua TRUE 1100.00
1144 Guatemala TRUE 1901.64
composite_drop_flag
TRUE
TRUE
TRUE
TRUE
TRUE
TRUE
TRUE
All corrected rows are altitude-flagged: TRUE
altitude_corrected_m is NA for all non-corrected rows: TRUE
Column classes:
row_id country_of_origin
"numeric" "character"
total_cup_points altitude_outlier_flag
"numeric" "logical"
altitude_corrected_m defect_outlier_flag
"numeric" "logical"
multivariate_grade_outlier_flag high_leverage_flag
"logical" "logical"
high_cook_flag composite_drop_flag
"logical" "logical"cd /app && echo "=== Deliverables in /app/outputs ===" && ls -1 outputs/ && echo "" && echo "=== analysis.R is self-contained (writes all 4) ===" && grep -cE 'write_csv|write_json|saveRDS|writeLines' analysis.R && echo "write calls found" && echo "" && echo "=== Headline numbers ===" && Rscript -e 'library(jsonlite); j<-fromJSON("outputs/outlier_report.json",simplifyVector=FALSE); cat(sprintf("n_input=%d -> n_after=%d\nalt fence [%.1f, %.1f]m | flagged=%d corrected=%d\ndefects flagged=%d/%d/%d\nMV: %d grades, thr=%.2f, flagged=%d\ninfluence: %s (AICs %.1f/%.1f/%.1f) lev=%d cook=%d both=%d\ncountry: %d eval, %d moved>=2\nsensitivity: flavor %.3f->%.3f (%.1f%%) CI[%.1f,%.1f] reps=%d\n", j$n_input,j$n_after_zero_score_drop, j$altitude$lower_bound_m,j$altitude$upper_bound_m,j$altitude$n_flagged,j$altitude$n_unit_corrected, j$defects$category_one_defects$n_flagged,j$defects$category_two_defects$n_flagged,j$defects$quakers$n_flagged, length(j$multivariate_grades$columns_used),j$multivariate_grades$threshold,j$multivariate_grades$n_flagged, j$influence$selected_model,j$influence$aic_grades_only,j$influence$aic_grades_altitude,j$influence$aic_grades_defects,j$influence$n_high_leverage,j$influence$n_high_cook,j$influence$n_both, j$country_ranking$n_countries_evaluated,j$country_ranking$n_rank_changes_ge_2, j$sensitivity$coef_flavor_full,j$sensitivity$coef_flavor_dropped,j$sensitivity$delta_pct,j$sensitivity$delta_pct_ci_lower,j$sensitivity$delta_pct_ci_upper,j$sensitivity$n_bootstrap_replicates))'=== Deliverables in /app/outputs === audit_memo.md coffee_ols_selected.rds outlier_flags.csv outlier_report.json === analysis.R is self-contained (writes all 4) === 4 write calls found === Headline numbers === n_input=1339 -> n_after=1338 alt fence [357.4, 4923.8]m | flagged=51 corrected=7 defects flagged=14/14/10 MV: 7 grades, thr=24.32, flagged=40 influence: grades_defects (AICs 4459.4/4461.3/4418.7) lev=110 cook=66 both=35 country: 21 eval, 2 moved>=2 sensitivity: flavor 2.085->1.653 (-20.7%) CI[-44.3,6.2] reps=500
No file changes, the agent left the workspace untouched.
# SyncVals normalized verifier summary # task: coffee-ratings-outliers # attempt: 2 status: completed raw_score: 1.0 # reward rule: binary: score 1.0 == 1.0 # => reward = 1.0 (PASS) # # 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). ---------------------------------------------------------------------- === Deliverables in /app/outputs === audit_memo.md coffee_ols_selected.rds outlier_flags.csv outlier_report.json === analysis.R is self-contained (writes all 4) === 4 write calls found === Headline numbers === n_input=1339 -> n_after=1338 alt fence [357.4, 4923.8]m | flagged=51 corrected=7 defects flagged=14/14/10 MV: 7 grades, thr=24.32, flagged=40 influence: grades_defects (AICs 4459.4/4461.3/4418.7) lev=110 cook=66 both=35 country: 21 eval, 2 moved>=2 sensitivity: flavor 2.085->1.653 (-20.7%) CI[-44.3,6.2] reps=500
Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_2d04af4db672466c. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_2d04af4db672466c · verifier authoritative; classifier explanatory.