coffee-ratings-outliers

cd /app && ls -la && echo "---DATA---" && ls -la data/ && echo "---OUTPUTS---" && ls -la outputs/ 2>/dev/null || echo "no outputs dir"

cd /app && cat analysis.R

cd /app && head -3 data/coffee_ratings.csv && echo "===COLS===" && head -1 data/coffee_ratings.csv | tr ',' '\n' | nl

cd /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>&1 | head -20

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")
g_all <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
# zero-score rows
z <- raw %>% filter(total_cup_points == 0)
cat("rows with total_cup_points==0:", nrow(z), "\n")
# all grades zero
allz <- raw %>% filter(if_all(all_of(g_all), ~ .x == 0))
cat("rows with all 10 grades == 0:", nrow(allz), "\n")
cat("total_cup_points of those:", allz$total_cup_points, "\n")
cat("summary total_cup_points:\n"); print(summary(raw$total_cup_points))
cat("how many total_cup_points < 50:", sum(raw$total_cup_points < 50, na.rm=TRUE), "\n")
print(raw %>% filter(total_cup_points < 50) %>% select(total_cup_points, aroma, flavor, aftertaste, acidity, body, balance, cupper_points))
' 2>&1 | head -40

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)
cat("=== altitude_mean_meters ===\n")
alt <- dat$altitude_mean_meters
cat("NA count:", sum(is.na(alt)), " n:", length(alt), "\n")
print(summary(alt))
cat("n positive:", sum(!is.na(alt) & alt>0), " n zero:", sum(!is.na(alt)& alt==0),"\n")
cat("Extreme high values (>8000):\n")
print(dat %>% filter(altitude_mean_meters>8000) %>% select(altitude, altitude_mean_meters) %>% head(20))
cat("Number >8000:", sum(alt>8000, na.rm=TRUE), "\n")
cat("=== log10 fence k=3 ===\n")
la <- log10(alt[!is.na(alt) & alt>0])
q <- quantile(la, c(0.25,0.75)); iqr <- q[2]-q[1]
lo <- q[1]-3*iqr; hi <- q[2]+3*iqr
cat("log lo/hi:", lo, hi, " meters:", 10^lo, 10^hi, "\n")
' 2>&1 | head -40

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
la <- log10(alt[!is.na(alt) & alt>0])
q <- quantile(la, c(0.25,0.75)); iqr <- q[2]-q[1]
lo <- 10^(q[1]-3*iqr); hi <- 10^(q[2]+3*iqr)
cat("fence meters:", lo, hi, "\n")
flag <- !is.na(alt) & (alt < lo | alt > hi)
cat("n flagged:", sum(flag), "\n")
fr <- dat %>% filter(flag) %>% select(altitude, altitude_mean_meters) %>% arrange(altitude_mean_meters)
print(fr, n=60)
' 2>&1 | head -70

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
la <- log10(alt[!is.na(alt) & alt>0])
q <- quantile(la, c(0.25,0.75)); iqr <- q[2]-q[1]
lo <- 10^(q[1]-3*iqr); hi <- 10^(q[2]+3*iqr)
flag <- !is.na(alt) & (alt < lo | alt > hi)
idx <- which(flag)
first_tok <- function(s){ m <- regmatches(s, regexpr("[0-9]+\\.?[0-9]*", s)); if(length(m)==0) NA_real_ else as.numeric(m) }
res <- lapply(idx, function(i){
  s <- dat$altitude[i]; tok <- first_tok(s)
  cand <- c(tok/10, tok/100, tok)
  inb <- cand[cand>=lo & cand<=hi & is.finite(cand)]
  corr <- if(length(inb)>0) inb[1] else NA_real_
  data.frame(altitude=s, mean_m=alt[i], tok=tok, corrected=corr)
})
res <- do.call(rbind, res)
print(res, row.names=FALSE)
cat("n_flagged:", sum(flag), " n_unit_corrected:", sum(!is.na(res$corrected)), "\n")
' 2>&1 | head -70

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]]
  cat("===",c,"=== NA:", sum(is.na(x)), "\n")
  cat("  pct zero:", round(mean(x==0,na.rm=TRUE),3), "\n")
  print(quantile(x, c(.5,.75,.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 IQR3 upper:", up, " n>up:", sum(x>up,na.rm=TRUE), "\n")
  cat("  table of top values:\n"); print(tail(sort(table(x)),8))
}
' 2>&1 | head -60

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
n <- nrow(dat)
for (c in c("category_one_defects","category_two_defects","quakers")){
  x <- dat[[c]]
  for (p in c(0.98,0.99)){
    thr <- as.numeric(quantile(x, p, na.rm=TRUE))
    nf <- sum(x>thr, na.rm=TRUE)
    cat(sprintf("%-22s p%.0f thr=%.3f  n>thr=%d (%.2f%%)\n", c, p*100, thr, nf, 100*nf/n))
  }
}
cat("\n--- combined OR across columns at p99 ---\n")
f <- rep(FALSE,n)
for (c in c("category_one_defects","category_two_defects","quakers")){
  x <- dat[[c]]; thr <- as.numeric(quantile(x,0.99,na.rm=TRUE)); f <- f | (!is.na(x)&x>thr)
}
cat("union flagged p99:", sum(f), sprintf("(%.2f%%)\n",100*sum(f)/n))
' 2>&1 | head -30

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")
clean <- c("uniformity","clean_cup","sweetness")
cat("=== 7 primary ===\n")
for(c in g7) cat(sprintf("%-14s sd=%.3f  min=%.2f\n",c,sd(dat[[c]]),min(dat[[c]])))
cat("=== cleanliness (near-constant) ===\n")
for(c in clean){ x<-dat[[c]]; cat(sprintf("%-12s sd=%.3f  pct==10: %.3f  min=%.2f\n",c,sd(x),mean(x==10),min(x))) }
# condition number of full-10 cov vs 7
G10<-as.matrix(dat[,c(g7,clean)]); G7<-as.matrix(dat[,g7])
cc<-complete.cases(G10)
cat("cond(cov 10):", kappa(cov(G10[cc,]),exact=TRUE), "\n")
cat("cond(cov 7):", kappa(cov(G7[complete.cases(G7),]),exact=TRUE), "\n")
# mahalanobis on 7
G<-G7[complete.cases(G7),]; md<-mahalanobis(G,colMeans(G),cov(G))
thr<-qchisq(0.999,df=7)
cat("thr df7 .999:",thr," n_flag:",sum(md>thr),"\n")
' 2>&1 | head -40

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")
cc <- complete.cases(dat[,g7])
cat("rows complete on 7 grades:", sum(cc), "of", nrow(dat), "\n")
fit_df <- dat[cc,]
# altitude corrected vector (reuse earlier logic)
alt <- dat$altitude_mean_meters
la <- log10(alt[!is.na(alt)&alt>0]); 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<lo|alt>hi)
first_tok<-function(s){m<-regmatches(s,regexpr("[0-9]+\\.?[0-9]*",s));if(length(m)==0)NA_real_ else as.numeric(m)}
corr<-rep(NA_real_,nrow(dat))
for(i in which(flag)){tok<-first_tok(dat$altitude[i]);cand<-c(tok/10,tok/100,tok);inb<-cand[is.finite(cand)&cand>=lo&cand<=hi];if(length(inb)>0)corr[i]<-inb[1]}
alt_use <- ifelse(!is.na(corr), corr, alt)
log_elev <- log10(alt_use)
m0 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points, data=fit_df)
# grades_altitude
fa <- fit_df; le <- log_elev[cc]; le[!is.finite(le)] <- mean(le[is.finite(le)]); fa$log_elev<-le
m1 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+log_elev, data=fa)
# grades_defects
fd <- fit_df
for(c in c("category_one_defects","category_two_defects","quakers")){x<-fd[[c]];x[is.na(x)]<-mean(x,na.rm=TRUE);fd[[paste0("l_",c)]]<-log1p(x)}
m2 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l_category_one_defects+l_category_two_defects+l_quakers, data=fd)
cat(sprintf("AIC grades_only=%.4f\nAIC grades_altitude=%.4f\nAIC grades_defects=%.4f\n", AIC(m0),AIC(m1),AIC(m2)))
cat("resid sd m0:", summary(m0)$sigma, "\n")
' 2>&1 | head -30

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")
fd <- dat
for(c in c("category_one_defects","category_two_defects","quakers")){x<-fd[[c]];x[is.na(x)]<-mean(x,na.rm=TRUE);fd[[paste0("l_",c)]]<-log1p(x)}
m2 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l_category_one_defects+l_category_two_defects+l_quakers, data=fd)
n<-nobs(m2); p<-length(coef(m2))
cat("n=",n," p=",p,"\n")
lev<-hatvalues(m2); ck<-cooks.distance(m2)
lev_thr<-2*p/n; ck_thr<-4/n
cat("lev_thr=",lev_thr," ck_thr=",ck_thr,"\n")
cat("n_high_lev=",sum(lev>lev_thr)," n_high_cook=",sum(ck>ck_thr)," n_both=",sum(lev>lev_thr & ck>ck_thr),"\n")
top<-order(-ck)[1:10]
cat("top10 cook row_ids (0-idx):", paste(top-1,collapse=","), "\n")
cat("top10 cook values:", paste(round(ck[top],3),collapse=","), "\n")
' 2>&1 | head -20

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)
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)
cat("n_countries:",nrow(ct)," n_changes_ge2:",sum(ct$d>=2)," max_change:",max(ct$d),"\n")
' 2>&1 | head -50

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")
# altitude flag
alt<-dat$altitude_mean_meters; la<-log10(alt[!is.na(alt)&alt>0]);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<lo|alt>hi)
# mv flag
G<-as.matrix(dat[,g7]); md<-mahalanobis(G,colMeans(G),cov(G)); flag_mv<-md>qchisq(0.999,df=7)
# influence (selected grades_defects)
fd<-dat; for(c in c("category_one_defects","category_two_defects","quakers")){x<-fd[[c]];x[is.na(x)]<-mean(x,na.rm=TRUE);fd[[paste0("l_",c)]]<-log1p(x)}
m2<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l_category_one_defects+l_category_two_defects+l_quakers,data=fd)
n<-nobs(m2);p<-length(coef(m2));lev<-hatvalues(m2);ck<-cooks.distance(m2)
hl<-lev>2*p/n; hc<-ck>4/n; infl<-hl&hc
composite<-flag_alt|flag_mv|infl
cat("n_alt:",sum(flag_alt)," n_mv:",sum(flag_mv)," n_infl(both):",sum(infl)," n_composite:",sum(composite),"\n")
f<-function(d)lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=d)
cf_full<-coef(f(dat))["flavor"]; cf_drop<-coef(f(dat[!composite,]))["flavor"]
dp<-100*(cf_drop-cf_full)/cf_full
cat(sprintf("coef_full=%.5f coef_drop=%.5f delta_pct=%.5f\n",cf_full,cf_drop,dp))
# bootstrap
set.seed(20260512); B<-500; nn<-nrow(dat); deltas<-numeric(0); succ<-0
for(b in 1:B){
  idx<-sample.int(nn,nn,replace=TRUE); db<-dat[idx,]; cb<-composite[idx]
  if(sum(!cb)<10) next
  ff<-tryCatch(f(db),error=function(e)NULL); fdr<-tryCatch(f(db[!cb,]),error=function(e)NULL)
  if(is.null(ff)||is.null(fdr)) next
  c1<-coef(ff)["flavor"]; c2<-coef(fdr)["flavor"]
  if(is.na(c1)||is.na(c2)||c1==0) next
  deltas<-c(deltas,100*(c2-c1)/c1); succ<-succ+1
}
ci<-quantile(deltas,c(.025,.975))
cat(sprintf("boot succ=%d ci=[%.4f, %.4f]\n",succ,ci[1],ci[2]))
' 2>&1 | head -20

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")
alt<-dat$altitude_mean_meters; la<-log10(alt[!is.na(alt)&alt>0]);q<-quantile(la,c(.25,.75));iqr<-q[2]-q[1]
lo<-10^(q[1]-3*iqr);hi<-10^(q[2]+3*iqr); A<-!is.na(alt)&(alt<lo|alt>hi)
G<-as.matrix(dat[,g7]); M<-mahalanobis(G,colMeans(G),cov(G))>qchisq(0.999,df=7)
D<-rep(FALSE,nrow(dat)); for(c in c("category_one_defects","category_two_defects","quakers")){x<-dat[[c]];thr<-as.numeric(quantile(x,0.99,na.rm=TRUE));D<-D|(!is.na(x)&x>thr)}
fd<-dat; for(c in c("category_one_defects","category_two_defects","quakers")){x<-fd[[c]];x[is.na(x)]<-mean(x,na.rm=TRUE);fd[[paste0("l_",c)]]<-log1p(x)}
m2<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l_category_one_defects+l_category_two_defects+l_quakers,data=fd)
n<-nobs(m2);p<-length(coef(m2));I<-(hatvalues(m2)>2*p/n)&(cooks.distance(m2)>4/n)
tot<-as.integer(A)+as.integer(M)+as.integer(D)+as.integer(I)
cat("by count 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("total flagged(any of 4):",sum(tot>0),"\n")
' 2>&1 | head

command: create
path: /app/analysis.R
file_text: #!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-aware rewrite
# ===================================================
# Each column is handled with a method that matches its distribution instead of
# the draft's one-size-fits-all meter-scale Tukey fences / all-column Mahalanobis
# / abs(rstandard)>2 influence / raw country means. See audit_memo.md for the
# per-block rationale. Re-running this script from an empty /app/outputs/
# regenerates every artifact.

suppressPackageStartupMessages({
  library(dplyr)
  library(readr)
  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")
CLEAN3  <- c("uniformity", "clean_cup", "sweetness")
DEFECTS <- c("category_one_defects", "category_two_defects", "quakers")

raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)

## ---- 1. Sentinel drop: the withdrawn lot with every grade == 0 -------------
all_grades <- c(GRADES7, CLEAN3)
zero_row <- rowSums(raw[, all_grades] == 0, na.rm = TRUE) == length(all_grades)
dat <- raw[!zero_row, , drop = FALSE]
n_after <- nrow(dat)
dat$row_id <- as.integer(seq_len(n_after) - 1L)
N <- n_after

## ---- 2. Altitude: Tukey fence (k=3) on log10 scale, unit correction --------
alt <- dat$altitude_mean_meters
la  <- log10(alt[!is.na(alt) & alt > 0])
qa  <- quantile(la, c(0.25, 0.75), names = FALSE)
iqr_a <- qa[2] - qa[1]
log_lo <- qa[1] - 3 * iqr_a
log_hi <- qa[2] + 3 * iqr_a
alt_lo_m <- 10 ^ log_lo
alt_hi_m <- 10 ^ log_hi
flag_alt <- !is.na(alt) & (alt < alt_lo_m | alt > alt_hi_m)
n_flag_alt <- sum(flag_alt)

first_numeric_token <- function(s) {
  if (is.na(s)) return(NA_real_)
  m <- regmatches(s, regexpr("[0-9]+\\.?[0-9]*", s))
  if (length(m) == 0 || m == "") NA_real_ else as.numeric(m)
}
alt_corrected <- rep(NA_real_, N)
for (i in which(flag_alt)) {
  tok <- first_numeric_token(dat$altitude[i])
  if (is.na(tok)) next
  cand <- c(tok / 10, tok / 100, tok)          # divide-by-10, divide-by-100, as-is
  inb  <- cand[is.finite(cand) & cand >= alt_lo_m & cand <= alt_hi_m]
  if (length(inb) > 0) alt_corrected[i] <- inb[1]
}
n_unit_corrected <- sum(!is.na(alt_corrected))

## ---- 3. Defects: upper-tail (99th pctile) thresholds, robust to zero mass ---
defect_flag <- function(x, p = 0.99) {
  thr <- as.numeric(quantile(x, p, na.rm = TRUE))
  list(upper_threshold = thr, flag = !is.na(x) & x > thr)
}
d_res <- lapply(DEFECTS, function(col) defect_flag(dat[[col]]))
names(d_res) <- DEFECTS
flag_defect <- Reduce(`|`, lapply(d_res, `[[`, "flag"))

## ---- 4. Multivariate grades: classical Mahalanobis on 7 primary attributes --
G <- as.matrix(dat[, GRADES7])
gc <- complete.cases(G)
mu <- colMeans(G[gc, , drop = FALSE])
S  <- cov(G[gc, , drop = FALSE])
md2 <- rep(NA_real_, N)
md2[gc] <- mahalanobis(G[gc, , drop = FALSE], mu, S)
mv_thresh <- qchisq(0.999, df = length(GRADES7))
flag_mv <- !is.na(md2) & md2 > mv_thresh
n_flag_mv <- sum(flag_mv)

## ---- 5. Influence: AIC model selection + leverage & Cook's D ---------------
fit_rows <- which(complete.cases(dat[, GRADES7]))
fdf <- dat[fit_rows, , drop = FALSE]

# elevation predictor: corrected metres where recovered, else raw mean metres
alt_use  <- ifelse(!is.na(alt_corrected), alt_corrected, alt)
log_elev_full <- log10(alt_use)
le <- log_elev_full[fit_rows]
le[!is.finite(le)] <- mean(le[is.finite(le)])
fdf$log_elev <- le

for (col in DEFECTS) {
  x <- fdf[[col]]
  x[is.na(x)] <- mean(x, na.rm = TRUE)
  fdf[[paste0("ld_", col)]] <- log1p(x)
}

f_grades <- reformulate(GRADES7, response = "total_cup_points")
f_alt    <- reformulate(c(GRADES7, "log_elev"), response = "total_cup_points")
f_def    <- reformulate(c(GRADES7, paste0("ld_", DEFECTS)), response = "total_cup_points")

fit_grades_only     <- lm(f_grades, data = fdf)
f
…[truncated]

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1 | tail -15 && echo "---FILES---" && ls -la outputs/

command: str_replace
path: /app/analysis.R
old_str: n_boot <- length(boot_deltas)
ci <- if (n_boot > 0) quantile(boot_deltas, c(0.025, 0.975), names = FALSE)
      else c(NA_real_, NA_real_)
new_str: n_boot <- length(boot_deltas)
ci <- if (n_boot > 0) {
  quantile(boot_deltas, c(0.025, 0.975), names = FALSE)
} else {
  c(NA_real_, NA_real_)
}

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1 | tail -15 && echo "---FILES---" && ls -la outputs/

cd /app && cat outputs/outlier_report.json

cd /app && 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))
cat("model class:",class(fit)," terms:",paste(names(coef(fit)),collapse=","),"\n")
cat("n=",n," p=",p,"\n")
lev<-hatvalues(fit); ck<-cooks.distance(fit)
hl<-lev>2*p/n; hc<-ck>4/n
cat("reload n_high_lev=",sum(hl)," n_high_cook=",sum(hc)," both=",sum(hl&hc),"\n")
# reconcile: fit rows == all 1338 rows in order
cat("CSV nrow=",nrow(csv)," match lev:",all(csv$high_leverage_flag==hl),
    " match cook:",all(csv$high_cook_flag==hc),"\n")
# composite reconcile
comp <- csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag | (csv$high_leverage_flag & csv$high_cook_flag)
cat("composite matches:",all(comp==csv$composite_drop_flag)," n_comp=",sum(comp),"\n")
# top10 cook reconcile
top10<-order(ck,decreasing=TRUE)[1:10]-1
cat("top10 from reload:",paste(top10,collapse=","),"\n")
cat("top10 from json  :",paste(unlist(js$influence$top10_cook_row_ids),collapse=","),"\n")
# altitude corrected only non-NA for flagged rows & in fence
lo<-js$altitude$lower_bound_m; hi<-js$altitude$upper_bound_m
cm<-csv$altitude_corrected_m; nn_c<-sum(!is.na(cm))
cat("n corrected:",nn_c," all in fence:",all(cm[!is.na(cm)]>=lo & cm[!is.na(cm)]<=hi),
    " all flagged:",all(csv$altitude_outlier_flag[!is.na(cm)]),"\n")
' 2>&1 | head -30

cd /app && Rscript -e '
suppressPackageStartupMessages({library(jsonlite); library(readr)})
js <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
req_top <- c("n_input","n_after_zero_score_drop","methods","altitude","defects",
  "multivariate_grades","influence","country_ranking","sensitivity","flag_interactions")
cat("top keys ok:", all(req_top %in% names(js)), "\n")
cat("methods keys:", paste(names(js$methods),collapse=","), "\n")
cat("influence keys:", paste(names(js$influence),collapse=","), "\n")
cat("sens keys:", paste(names(js$sensitivity),collapse=","), "\n")
fi <- js$flag_interactions
cat("flag_interactions keys (",length(fi),"):", paste(names(fi),collapse=","), "\n")
cat("per_country len:", length(js$country_ranking$per_country), "\n")
cat("cols_used len:", length(js$multivariate_grades$columns_used), "\n")
# types
cat("n_input is int-like:", js$n_input==as.integer(js$n_input), "\n")
cat("aic types numeric:", is.numeric(js$influence$aic_grades_defects), "\n")
# CSV header check
hdr <- names(read_csv("outputs/outlier_flags.csv", show_col_types=FALSE, n_max=1))
exp <- c("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")
cat("CSV header exact:", identical(hdr, exp), "\n")
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
cat("defect union in CSV:", sum(csv$defect_outlier_flag), "\n")
cat("row_id 0-indexed contiguous:", identical(csv$row_id, 0:(nrow(csv)-1)), "\n")
' 2>&1 | head -30

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
cat("class row_id:", class(csv$row_id), "\n")
cat("values equal 0..n-1:", all(csv$row_id == 0:(nrow(csv)-1)), "\n")
cat("first/last:", head(csv$row_id,3), "...", tail(csv$row_id,3), "\n")
# check post-drop input order preserved: total_cup_points matches raw minus sentinel
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
ag <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points","uniformity","clean_cup","sweetness")
keep <- rowSums(raw[,ag]==0, na.rm=TRUE) != length(ag)
cat("order preserved:", all(csv$total_cup_points == raw$total_cup_points[keep]), "\n")
' 2>&1 | head

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
r <- raw[raw$altitude=="518" & !is.na(raw$altitude), c("altitude","unit_of_measurement","altitude_mean_meters")]
print(unique(r)); cat("518*0.3048 =", 518*0.3048, "(feet->m)\n")
' 2>&1 | head
echo "=== determinism: re-run and diff ==="
cd /app && cp outputs/outlier_report.json /tmp/run1.json && Rscript analysis.R >/dev/null 2>&1 && diff -q /tmp/run1.json outputs/outlier_report.json && echo "JSON identical across runs" || echo "DIFFERS"

command: create
path: /app/outputs/audit_memo.md
file_text: # Coffee Quality Outlier Audit , method memo

The draft applied the same off-the-shelf summary to every column. Each section
below says why that fails for the column's actual shape and what replaced it.
All counts are from the re-run (`n_input = 1339`, `n_after = 1338`).

## Data prep
The draft dropped `total_cup_points > 0`, which happens to remove the right row
but encodes the wrong rule. The documented sentinel is the withdrawn lot whose
**ten grades are all zero**; I drop on that condition (exactly one row, the same
one) so the intent is explicit. That lot is removed *before* any quantile, mean
or covariance is computed , left in, its zeros would drag every fence, the grade
covariance and all country means. Surviving lots get a 0-indexed `row_id` in
post-drop input order so the CSV, the model rows and the JSON all align.

## Altitude
`altitude_mean_meters` is heavily right-skewed (median 1311 m, max 190164 m), so
the draft's Tukey fence on the **raw metre scale** is stretched by the long tail
and is useless at the low end; it also never attempts unit recovery. I take
`log10` of the positive values and build the Tukey fence with `k = 3`, then
back-transform: **[357.4 m, 4923.8 m]**, flagging **51** rows. For each flagged
row I test power-of-ten corrections on the first numeric token of the raw
`altitude` string , `÷10`, then `÷100`, then as-is , keeping the first candidate
inside the fence. This recovers **7** lots: e.g. `11000 metros → 1100`,
`190164 → 1901.64`, and three `518` lots stored as feet (`518 ft` mis-parsed to
158 m) recovered as-is to 518 m. Rows with no in-bounds candidate (`"1"`,
`"200 m"`) keep `altitude_corrected_m = NA`.

## Defects
The draft's Tukey IQR×3 on raw counts is the worst possible choice here. These
columns are mass-at-zero (85% / 28% / 93% zeros), so for `category_one_defects`
and `quakers` **Q1 = Q3 = 0**: the fence collapses to 0 and flags *every* nonzero
value , 202 and 94 rows (≈15% and ≈7%), the opposite of "genuinely extreme." Any
IQR/MAD fence degenerates once a majority ties at zero. I instead cut each column
at its **99th percentile** (7.63 / 26.63 / 5), flagging **14 / 14 / 10** lots
(~1% each) , matching the stated sparse top 1–2% tail. A row is a defect outlier
if any column trips (union = 35).

## Multivariate grades
The draft ran Mahalanobis on **all ten** grade columns, including
`uniformity`, `clean_cup`, `sweetness` , near-constants pinned at exactly 10 for
86–91% of lots. Those degenerate directions inflate the covariance condition
number and let microscopic wobble in the pinned columns dominate the distance, so
the statistic stops measuring joint *flavour* structure. I use the **seven
primary SCA attributes** with the classical sample mean/covariance (not MCD, per
spec) and a chi-square cutoff at 99.9% with df = 7 (**threshold 24.32**),
flagging **40** lots.

## Influence
`abs(rstandard) > 2` measures large *residuals* (outliers in `y`), not
*influence*, and hard-codes one model. I fit the three specified candidates on
the 1338 grade-complete rows and select the lowest finite **AIC**:
`grades_defects` (4418.7) < `grades_only` (4459.4) < `grades_altitude` (4461.3).
`total_cup_points` embeds cleanliness deductions, which the `log1p` defect terms
proxy , hence the win. On that fit I flag high leverage (`> 2p/n = 0.0164`) and
high Cook's D (`> 4/n = 0.00299`) **separately** (110 and 66 lots); the influence
flag requires **both** (35). The fitted model is saved to
`coffee_ols_selected.rds` for re-verification.

## Country ranking
The draft reported only raw means (trimmed mean set equal to raw), so there was
no robustness check at all. For the **21** countries with ≥10 lots I compare the
raw mean against a **10% trimmed mean**. The leaderboard is mostly stable: only
**2** countries move ≥2 positions , Taiwan (14→16) and Honduras (20→18, where
trimming removes a cluster of low-scoring lots) , with a maximum shift of 2. The
…[truncated]

cd /app && echo "word count: $(wc -w < outputs/audit_memo.md)" && echo "--- headings ---" && grep '^## ' outputs/audit_memo.md && echo "--- required check ---" && for h in "Data prep" "Altitude" "Defects" "Multivariate grades" "Influence" "Country ranking" "Sensitivity" "Composite flags"; do grep -q "^## $h\$" outputs/audit_memo.md && echo "OK: $h" || echo "MISSING: $h"; done

command: str_replace
path: /app/analysis.R
old_str: 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), n_boot))
new_str: ## ---- Audit memo (regenerated on every run) --------------------------------
memo <- r"=====(# Coffee Quality Outlier Audit -- method memo

The draft applied the same off-the-shelf summary to every column. Each section
below says why that fails for the column's actual shape and what replaced it.
All counts are from the re-run (`n_input = 1339`, `n_after = 1338`).

## Data prep
The draft dropped `total_cup_points > 0`, which happens to remove the right row
but encodes the wrong rule. The documented sentinel is the withdrawn lot whose
**ten grades are all zero**; I drop on that condition (exactly one row, the same
one) so the intent is explicit. That lot is removed *before* any quantile, mean
or covariance is computed -- left in, its zeros would drag every fence, the grade
covariance and all country means. Surviving lots get a 0-indexed `row_id` in
post-drop input order so the CSV, the model rows and the JSON all align.

## Altitude
`altitude_mean_meters` is heavily right-skewed (median 1311 m, max 190164 m), so
the draft's Tukey fence on the **raw metre scale** is stretched by the long tail
and is useless at the low end; it also never attempts unit recovery. I take
`log10` of the positive values and build the Tukey fence with `k = 3`, then
back-transform: **[357.4 m, 4923.8 m]**, flagging **51** rows. For each flagged
row I test power-of-ten corrections on the first numeric token of the raw
`altitude` string -- /10, then /100, then as-is -- keeping the first candidate
inside the fence. This recovers **7** lots: e.g. `11000 metros -> 1100`,
`190164 -> 1901.64`, and three `518` lots stored as feet (`518 ft` mis-parsed to
158 m) recovered as-is to 518 m. Rows with no in-bounds candidate (`"1"`,
`"200 m"`) keep `altitude_corrected_m = NA`.

## Defects
The draft's Tukey IQR*3 on raw counts is the worst possible choice here. These
columns are mass-at-zero (85% / 28% / 93% zeros), so for `category_one_defects`
and `quakers` **Q1 = Q3 = 0**: the fence collapses to 0 and flags *every* nonzero
value -- 202 and 94 rows (~15% and ~7%), the opposite of "genuinely extreme." Any
IQR/MAD fence degenerates once a majority ties at zero. I instead cut each column
at its **99th percentile** (7.63 / 26.63 / 5), flagging **14 / 14 / 10** lots
(~1% each) -- matching the stated sparse top 1-2% tail. A row is a defect outlier
if any column trips (union = 35).

## Multivariate grades
The draft ran Mahalanobis on **all ten** grade columns, including
`uniformity`, `clean_cup`, `sweetness` -- near-constants pinned at exactly 10 for
86-91% of lots. Those degenerate directions inflate the covariance condition
number and let microscopic wobble in the pinned columns dominate the distance, so
the statistic stops measuring joint *flavour* structure. I use the **seven
primary SCA attributes** with the classical sample mean/covariance (not MCD, per
spec) and a chi-square cutoff at 99.9% with df = 7 (**threshold 24.32**),
flagging **40** lots.

## Influence
`abs(rstandard) > 2` measures large *residuals* (outliers in `y`), not
*influence*, and hard-codes one model. I fit the three specified candidates on
the 1338 grade-complete rows and select the lowest finite **AIC**:
`grades_defects` (4418.7) < `grades_only` (4459.4) < `grades_altitude` (4461.3).
`total_cup_points` embeds cleanliness deductions, which the `log1p` defect terms
proxy -- hence the win. On that fit I flag high leverage (`> 2p/n = 0.0164`) and
high Cook's D (`> 4/n = 0.00299`) **separately** (110 and 66 lots); the influence
flag requires **both** (35). The fitted model is saved to
`coffee_ols_selected.rds` for re-verification.

## Country ranking
The draft reported only raw means (trimmed mean set equal to raw), so there was
no robustness check at all. For the **21** countries with 
…[truncated]

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1 | tail -5 && echo "---FILES---" && ls -1 outputs/ && echo "---MEMO CHECK---" && echo "words: $(wc -w < outputs/audit_memo.md)" && grep -c '^## ' outputs/audit_memo.md && head -1 outputs/audit_memo.md

cd /app && echo "occurrences of 5+ dashes in memo content:" && grep -c -- '-----' analysis.R; echo "occurrences of )----- :" && grep -c ')-----' analysis.R

cd /app && echo ")--- count:" && grep -c ')---' analysis.R && echo "where 5-dashes occur:" && grep -n -- '-----' analysis.R | head

cd /app && printf 'lines with 5 dashes: '; grep -c -- '-----' analysis.R; printf 'lines with 4 dashes: '; grep -c -- '----' analysis.R; printf 'lines with )-----": '; grep -c ')-----"' analysis.R

cd /app && grep -n -- '-----' analysis.R

cd /app && Rscript -e 'x <- r"-----(a "quoted" and apostrophe'"'"'s test with -- dashes)-----"; cat(x, "\n")' 2>&1

cd /app && sed -i 's/memo <- r"=====(/memo <- r"-----(/; s/^)====="$/)-----"/' analysis.R && echo "--- delimiter lines now ---" && grep -n 'memo <- r"' analysis.R && grep -n '^)-----"$' analysis.R

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1 | tail -5 && echo "---FILES---" && ls -1 outputs/ && echo "---MEMO---" && echo "words: $(wc -w < outputs/audit_memo.md), headings: $(grep -c '^## ' outputs/audit_memo.md)" && head -1 outputs/audit_memo.md

cd /app && Rscript -e '
suppressPackageStartupMessages({library(jsonlite); library(readr)})
stopifnot(all(file.exists(file.path("outputs", c("outlier_report.json","coffee_ols_selected.rds","outlier_flags.csv","audit_memo.md")))))
js <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
fit <- readRDS("outputs/coffee_ols_selected.rds")
n<-nobs(fit); p<-length(coef(fit)); lev<-hatvalues(fit); ck<-cooks.distance(fit)
# 1 reconcile flags from reloaded model
ok_lev <- all((lev>2*p/n)==csv$high_leverage_flag)
ok_ck  <- all((ck>4/n)==csv$high_cook_flag)
ok_comp<- all((csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag))==csv$composite_drop_flag)
ok_top <- all(unlist(js$influence$top10_cook_row_ids)==(order(ck,decreasing=TRUE)[1:10]-1))
# 2 JSON internal consistency
ok_cnt <- js$influence$n_both==sum(csv$high_leverage_flag&csv$high_cook_flag)
ok_mv  <- js$multivariate_grades$n_flagged==sum(csv$multivariate_grade_outlier_flag)
ok_alt <- js$altitude$n_flagged==sum(csv$altitude_outlier_flag)
ok_def <- sum(sapply(js$defects,function(d)d$n_flagged))>=js$defects$quakers$n_flagged # sanity
ok_corr<- js$altitude$n_unit_corrected==sum(!is.na(csv$altitude_corrected_m))
# 3 country ranks consistent with reported means
pc <- js$country_ranking$per_country
rm_<-sapply(pc,function(x)x$raw_mean); tm_<-sapply(pc,function(x)x$trimmed_mean)
rr<-sapply(pc,function(x)x$raw_rank); tr<-sapply(pc,function(x)x$trimmed_rank)
ok_rank <- all(rr==rank(-rm_,ties.method="min")) && all(tr==rank(-tm_,ties.method="min"))
ok_chg  <- js$country_ranking$n_rank_changes_ge_2==sum(abs(rr-tr)>=2)
cat(sprintf("reload lev/cook: %s/%s  composite: %s  top10: %s\n",ok_lev,ok_ck,ok_comp,ok_top))
cat(sprintf("json counts mv/alt/corr/both: %s/%s/%s/%s\n",ok_mv,ok_alt,ok_corr,ok_cnt))
cat(sprintf("country ranks/changes: %s/%s\n",ok_rank,ok_chg))
cat(sprintf("selected=%s n=%d p=%d  flavor full/drop=%.4f/%.4f\n",js$influence$selected_model,n,p,js$sensitivity$coef_flavor_full,js$sensitivity$coef_flavor_dropped))
cat("ALL TRUE:", all(ok_lev,ok_ck,ok_comp,ok_top,ok_mv,ok_alt,ok_corr,ok_cnt,ok_rank,ok_chg),"\n")
' 2>&1 | tail -12

cd /app && python3 -c "
import json
d=json.load(open('outputs/outlier_report.json'))
print('JSON valid. top-level keys:', len(d))
print('n_input',d['n_input'],'n_after',d['n_after_zero_score_drop'])
print('alt bounds', round(d['altitude']['lower_bound_m'],2), round(d['altitude']['upper_bound_m'],2), 'flagged',d['altitude']['n_flagged'],'corrected',d['altitude']['n_unit_corrected'])
print('mv cols', d['multivariate_grades']['columns_used'])
print('aic', d['influence']['aic_grades_only'],d['influence']['aic_grades_altitude'],d['influence']['aic_grades_defects'])
print('sens', d['sensitivity'])
print('flag_inter sum 0-4 =', sum(d['flag_interactions'][k] for k in ['n_with_zero_flags','n_with_one_flag','n_with_two_flags','n_with_three_flags','n_with_four_flags']))
"
echo "--- CSV rows with corrected altitude ---"
awk -F, 'NR==1{for(i=1;i<=NF;i++)h[i]=\$i} NR>1 && \$5!=\"\"{print \$1,\$5}' outputs/outlier_flags.csv

cd /app && Rscript -e '
suppressPackageStartupMessages(library(readr))
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
sub <- csv[!is.na(csv$altitude_corrected_m), c("row_id","country_of_origin","altitude_outlier_flag","altitude_corrected_m")]
print(as.data.frame(sub), row.names=FALSE)
cat("all corrected rows are altitude-flagged:", all(sub$altitude_outlier_flag), "\n")
cat("NA corrected count:", sum(is.na(csv$altitude_corrected_m)), " non-NA:", sum(!is.na(csv$altitude_corrected_m)), "\n")
# show raw NA rendering in file
cat("--- sample raw lines (corrected vs not) ---\n")
' 2>&1
echo "--- grep a corrected row + a normal row in raw csv ---"
head -1 outputs/outlier_flags.csv
awk -F, 'NR>1 && $5!=""' outputs/outlier_flags.csv | head -3
echo "(normal row:)"; sed -n '3p' outputs/outlier_flags.csv

cd /app && rm -rf outputs && mkdir -p outputs && Rscript -e 'options(warn=1); source("analysis.R")' 2>&1 | grep -iE "warn|error|deficient|singular|NaN|collinear" || echo "No warnings/errors emitted."

cd /app && sed -n '1,120p' analysis.R

cd /app && sed -n '120,205p' analysis.R

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R && echo "=== OUTPUT FILES ===" && ls -la outputs/ && echo && echo "=== KEY RESULTS ===" && python3 -c "
import json
d=json.load(open('outputs/outlier_report.json'))
print(f\"n_input={d['n_input']}  n_after={d['n_after_zero_score_drop']}\")
a=d['altitude']; print(f\"altitude: fence=[{a['lower_bound_m']:.1f},{a['upper_bound_m']:.1f}] flagged={a['n_flagged']} corrected={a['n_unit_corrected']}\")
de=d['defects']; print('defects:', {k:(de[k]['upper_threshold'],de[k]['n_flagged']) for k in de})
m=d['multivariate_grades']; print(f\"mv: df7 thr={m['threshold']:.2f} flagged={m['n_flagged']} cols={len(m['columns_used'])}\")
i=d['influence']; print(f\"influence: {i['selected_model']} AIC={i['aic_grades_only']:.1f}/{i['aic_grades_altitude']:.1f}/{i['aic_grades_defects']:.1f} lev={i['n_high_leverage']} cook={i['n_high_cook']} both={i['n_both']}\")
c=d['country_ranking']; print(f\"country: eval={c['n_countries_evaluated']} changes>=2={c['n_rank_changes_ge_2']} max={c['max_abs_rank_change']}\")
s=d['sensitivity']; print(f\"sens: flavor {s['coef_flavor_full']:.3f}->{s['coef_flavor_dropped']:.3f} delta={s['delta_pct']:.1f}% CI=[{s['delta_pct_ci_lower']:.1f},{s['delta_pct_ci_upper']:.1f}] B={s['n_bootstrap_replicates']}\")
f=d['flag_interactions']; print('interactions 0-4:', [f[k] for k in ['n_with_zero_flags','n_with_one_flag','n_with_two_flags','n_with_three_flags','n_with_four_flags']])
"