coffee-ratings-outliers

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

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

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\n")
grades <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
# Find zero-score rows
allzero <- rowSums(raw[,grades] == 0, na.rm=TRUE) == length(grades)
cat("rows with ALL 10 grades == 0:", sum(allzero, na.rm=TRUE), "\n")
cat("rows with total_cup_points == 0:", sum(raw$total_cup_points == 0, na.rm=TRUE), "\n")
cat("rows with total_cup_points <= 0:", sum(raw$total_cup_points <= 0, na.rm=TRUE), "\n\n")
cat("Summary of grade columns:\n")
print(sapply(raw[,grades], function(x) c(min=min(x,na.rm=T), med=median(x,na.rm=T), max=max(x,na.rm=T), sd=round(sd(x,na.rm=T),3))))
'

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: NA count =", sum(is.na(alt)), " of", nrow(dat), "\n")
cat("n positive =", sum(!is.na(alt) & alt>0), "  n zero =", sum(!is.na(alt) & alt==0), "\n")
cat("quantiles (raw m):\n"); print(quantile(alt, c(0,.01,.25,.5,.75,.9,.99,1), na.rm=TRUE))
la <- log10(alt[!is.na(alt) & alt>0])
q <- quantile(la, c(.25,.75)); iqr <- q[2]-q[1]
lo <- q[1]-3*iqr; hi <- q[2]+3*iqr
cat("\nlog10 Q1,Q3 =", round(q,4), " IQR =", round(iqr,4), "\n")
cat("log10 fence = [", round(lo,4), ",", round(hi,4), "]\n")
cat("meters fence = [", round(10^lo,2), ",", round(10^hi,2), "]\n")
flg <- !is.na(alt) & (alt < 10^lo | alt > 10^hi)
cat("n flagged outside log fence =", sum(flg), "\n")
cat("\nFlagged altitude values (mean_m) and raw altitude string:\n")
fd <- dat[flg, c("altitude","altitude_mean_meters","unit_of_measurement")]
print(as.data.frame(fd), row.names=FALSE)
'

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(.25,.75)); iqr <- q[[2]]-q[[1]]
lo_m <- 10^(q[[1]]-3*iqr); hi_m <- 10^(q[[2]]+3*iqr)
cat("fence meters:", round(lo_m,4), round(hi_m,4), "\n")
flg <- !is.na(alt) & (alt < lo_m | alt > hi_m)
idx <- which(flg)
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 =", sum(flg), "  n_unit_corrected =", ncorr, "\n\n")
show <- data.frame(altitude=dat$altitude[idx], mean_m=alt[idx], corrected=corrected[idx])
print(show, row.names=FALSE)
'

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 (col in c("category_one_defects","category_two_defects","quakers")){
  x <- dat[[col]]
  cat("==", col, "== NA:", sum(is.na(x)), "\n")
  cat("  pct zero:", round(mean(x==0, na.rm=TRUE)*100,1), "%\n")
  print(table(x, useNA="ifany"))
  cat("  q50,75,90,95,98,99:", round(quantile(x, c(.5,.75,.9,.95,.98,.99), na.rm=TRUE),2), "\n")
  # Tukey raw
  q <- quantile(x, c(.25,.75), na.rm=TRUE); up <- q[[2]]+3*(q[[2]]-q[[1]])
  cat("  raw Tukey IQR3 upper:", up, " n>up:", sum(x>up, na.rm=TRUE), "\n\n")
}
'

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 (col in c("category_one_defects","category_two_defects","quakers")){
  x <- dat[[col]]; xn <- x[!is.na(x)]
  pos <- xn[xn>0]
  # Method A: Tukey k=3 on positive-only
  qa <- quantile(pos, c(.25,.75)); upA <- qa[[2]]+3*(qa[[2]]-qa[[1]])
  # Method B: Tukey k=1.5 on positive-only
  upB <- qa[[2]]+1.5*(qa[[2]]-qa[[1]])
  # Method C: 98th pctile of all
  upC <- quantile(xn, .98)
  cat(sprintf("%-22s n=%d posN=%d | posTukey3 up=%.1f nFlag=%d (%.1f%%) | posTukey1.5 up=%.1f nFlag=%d (%.1f%%) | p98=%.1f nFlag=%d (%.1f%%)\n",
    col, length(xn), length(pos),
    upA, sum(xn>upA), 100*sum(xn>upA)/n,
    upB, sum(xn>upB), 100*sum(xn>upB)/n,
    upC, sum(xn>upC), 100*sum(xn>upC)/n))
}
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G <- as.matrix(dat[,g7]); cc <- complete.cases(G)
cat("complete cases on 7 grades:", sum(cc), "of", nrow(dat), "\n")
mu <- colMeans(G[cc,]); S <- cov(G[cc,])
md2 <- mahalanobis(G[cc,], mu, S)
thr <- qchisq(0.999, df=7)
cat("threshold chisq(0.999, df=7) =", round(thr,4), "\n")
cat("n flagged (7-grade classical MD) =", sum(md2>thr), "\n")
# compare to all-10 (draft) 
g10 <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
G10 <- as.matrix(dat[,g10]); cc10 <- complete.cases(G10)
S10 <- cov(G10[cc10,])
cat("\ncondition number cov(7) =", round(kappa(S),1), " | cov(10) =", round(kappa(S10),1), "\n")
cat("sd of cleanliness cols:", round(sapply(dat[c("uniformity","clean_cup","sweetness")], sd, na.rm=TRUE),3), "\n")
'

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(as.matrix(dat[,g7]))
cat("rows complete on 7 grades:", sum(cc), "\n")
fit_df <- dat[cc,]
# altitude: corrected where available else mean; compute log10
alt <- fit_df$altitude_mean_meters
# approximate corrected: reuse fence correction quickly (use mean only for AIC test here)
logalt <- log10(alt); logalt[!is.finite(logalt)] <- NA
logalt[is.na(logalt)] <- mean(logalt, na.rm=TRUE)
fit_df$log_alt <- logalt
for (col in c("category_one_defects","category_two_defects","quakers")){
  x <- fit_df[[col]]; x[is.na(x)] <- mean(x, na.rm=TRUE); fit_df[[paste0("l1p_",col)]] <- log1p(x)
}
m1 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points, data=fit_df)
m2 <- update(m1, . ~ . + log_alt)
m3 <- lm(total_cup_points ~ aroma+flavor+aftertaste+acidity+body+balance+cupper_points+l1p_category_one_defects+l1p_category_two_defects+l1p_quakers, data=fit_df)
cat(sprintf("AIC grades_only=%.4f  grades_altitude=%.4f  grades_defects=%.4f\n", AIC(m1),AIC(m2),AIC(m3)))
cat("resid df / n used:", nobs(m1), nobs(m2), nobs(m3), "\n")
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types = FALSE)
dat <- raw %>% filter(total_cup_points > 0)
for (tr in c(0.1, 0.2)){
  ct <- dat %>% filter(!is.na(country_of_origin)) %>%
    group_by(country=country_of_origin) %>%
    summarise(n=n(), raw_mean=mean(total_cup_points), trimmed_mean=mean(total_cup_points, trim=tr), .groups="drop") %>%
    filter(n>=10) %>%
    mutate(raw_rank=rank(-raw_mean, ties.method="min"),
           trimmed_rank=rank(-trimmed_mean, ties.method="min"),
           chg=abs(raw_rank-trimmed_rank))
  cat(sprintf("trim=%.2f: n_countries=%d  n_changes_ge2=%d  max_chg=%d\n", tr, nrow(ct), sum(ct$chg>=2), max(ct$chg)))
}
cat("\nCountries with >=10 lots:", nrow(dat %>% filter(!is.na(country_of_origin)) %>% count(country_of_origin) %>% filter(n>=10)), "\n")
'

command: create
path: /app/analysis.R
file_text: #!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-aware rewrite
# ===================================================
# Each column is audited with a method that matches its distribution:
#   * altitude  -> Tukey fence on log10 scale (right-skewed, heavy tail)
#   * defects   -> 98th-percentile upper-tail cut (mass-at-zero counts)
#   * grades    -> classical Mahalanobis on the 7 primary SCA attributes only
#   * influence -> AIC-selected OLS, leverage>2p/n AND Cook's D>4/n
#   * country   -> raw mean vs 20% trimmed mean rank comparison
# The script is self-contained and reproduces every artifact from a clean
# /app/outputs/ directory.

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(readr)
  library(stringr)
  library(jsonlite)
  library(purrr)
})

DATA_PATH  <- "/app/data/coffee_ratings.csv"
OUTPUT_DIR <- "/app/outputs"
dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE)

GRADES7 <- c("aroma", "flavor", "aftertaste", "acidity",
             "body", "balance", "cupper_points")
GRADES10 <- c(GRADES7, "uniformity", "clean_cup", "sweetness")
DEFECTS <- c("category_one_defects", "category_two_defects", "quakers")

# ----------------------------------------------------------------------------
# 1. Load + sentinel drop
# ----------------------------------------------------------------------------
raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)

# The withdrawn submission has every one of the ten grades recorded as zero.
all_grades_zero <- rowSums(raw[, GRADES10] == 0, na.rm = TRUE) == length(GRADES10)
dat <- raw[!all_grades_zero, , drop = FALSE]
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L))

# ----------------------------------------------------------------------------
# 2. Altitude: Tukey fence on the log10 scale (k = 3), back-transformed
# ----------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
log_alt_pos <- log10(alt[!is.na(alt) & alt > 0])
qa <- quantile(log_alt_pos, c(0.25, 0.75), names = FALSE)
iqr_a <- qa[2] - qa[1]
lo_log <- qa[1] - 3 * iqr_a
hi_log <- qa[2] + 3 * iqr_a
lo_m <- 10^lo_log
hi_m <- 10^hi_log
flag_alt <- !is.na(alt) & (alt < lo_m | alt > hi_m)
n_flag_alt <- sum(flag_alt)

# Unit-slip recovery: on each flagged row, test power-of-ten corrections of the
# first numeric token in the raw `altitude` string (/10, then /100, then as-is)
# and keep the first candidate that lands inside the fence.
first_numeric_token <- function(s) suppressWarnings(as.numeric(str_extract(s, "[0-9]+\\.?[0-9]*")))
alt_corrected <- rep(NA_real_, nrow(dat))
for (i in which(flag_alt)) {
  tok <- first_numeric_token(dat$altitude[i])
  if (is.na(tok)) next
  candidates <- c(tok / 10, tok / 100, tok)
  in_fence <- candidates[candidates >= lo_m & candidates <= hi_m]
  if (length(in_fence) > 0) alt_corrected[i] <- in_fence[1]
}
n_unit_corrected <- sum(!is.na(alt_corrected))

# ----------------------------------------------------------------------------
# 3. Defects: 98th-percentile upper-tail cut per column (mass-at-zero counts)
# ----------------------------------------------------------------------------
defect_tail <- function(x) {
  thr <- quantile(x, 0.98, na.rm = TRUE, names = FALSE)
  flagged <- !is.na(x) & x > thr
  list(upper_threshold = thr, n = sum(flagged), flag_vec = flagged)
}
d_list <- lapply(DEFECTS, function(col) defect_tail(dat[[col]]))
names(d_list) <- DEFECTS
flag_defect <- Reduce(`|`, lapply(d_list, `[[`, "flag_vec"))

# ----------------------------------------------------------------------------
# 4. Multivariate grades: classical Mahalanobis on the 7 primary attributes
# ----------------------------------------------------------------------------
G <- as.matrix(dat[, GRADES7])
G_complete <- complete.cases(G)
mu <- colMeans(G[G_complete, ])
S  <- cov(G[G_complete, ])
md2 <- rep(NA_real_, nrow(dat))
md2[G_complete
…[truncated]

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1

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")
n <- nobs(fit); p <- length(coef(fit))
cat("RELOAD fit: class=", class(fit), " n=", n, " p=", p, " terms:", paste(names(coef(fit)), collapse=","), "\n")
lev <- hatvalues(fit); ck <- cooks.distance(fit)
lev_flag <- lev > 2*p/n; cook_flag <- ck > 4/n
# fit rows align 1:1 with csv rows (all complete)
cat("CSV rows =", nrow(csv), " fit rows =", length(lev), "\n")
cat("leverage reconcile:", all(lev_flag == csv$high_leverage_flag), " sum:", sum(lev_flag), "vs", sum(csv$high_leverage_flag), "\n")
cat("cook reconcile:", all(cook_flag == csv$high_cook_flag), " sum:", sum(cook_flag), "vs", sum(csv$high_cook_flag), "\n")
# composite reconcile
comp <- csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag | (csv$high_leverage_flag & csv$high_cook_flag)
cat("composite reconcile:", all(comp == csv$composite_drop_flag), "\n")
# top10 cook
ord <- order(ck, decreasing=TRUE)[1:10]
cat("top10 cook row_ids (recomputed):", paste(csv$row_id[ord], collapse=","), "\n")
cat("top10 in json:                   ", paste(js$influence$top10_cook_row_ids, collapse=","), "\n")
cat("top10 match:", all(csv$row_id[ord] == js$influence$top10_cook_row_ids), "\n")
# thresholds match json
cat("lev_thresh json=", js$influence$leverage_threshold, " recomputed=", 2*p/n, "\n")
cat("cook_thresh json=", js$influence$cook_threshold, " recomputed=", 4/n, "\n")
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite); library(dplyr)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types = FALSE)
js  <- fromJSON("outputs/outlier_report.json")
cat("CSV columns:\n"); print(names(csv))
cat("\nrow_id range:", min(csv$row_id), "-", max(csv$row_id), " unique:", length(unique(csv$row_id)), "\n")
cat("altitude_corrected_m non-NA count:", sum(!is.na(csv$altitude_corrected_m)), " (all within corrected-only?)\n")
cat("  corrected rows all flagged as altitude_outlier:", all(csv$altitude_outlier_flag[!is.na(csv$altitude_corrected_m)]), "\n")
cat("  corrected values:", paste(round(na.omit(csv$altitude_corrected_m),2), collapse=", "), "\n")
cat("  all corrected within fence [", round(js$altitude$lower_bound_m,2),",",round(js$altitude$upper_bound_m,2),"]:",
    all(na.omit(csv$altitude_corrected_m) >= js$altitude$lower_bound_m & na.omit(csv$altitude_corrected_m) <= js$altitude$upper_bound_m), "\n\n")
# defect reconcile using raw data thresholds
raw <- read_csv("/app/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")){
  thr <- js$defects[[col]]$upper_threshold
  nf <- sum(dat[[col]] > thr, na.rm=TRUE)
  cat(sprintf("%-22s thr=%-4g json_n=%d recompute(x>thr)=%d match=%s\n", col, thr, js$defects[[col]]$n_flagged, nf, js$defects[[col]]$n_flagged==nf))
}
# union defect
du <- (dat$category_one_defects > js$defects$category_one_defects$upper_threshold) |
      (dat$category_two_defects > js$defects$category_two_defects$upper_threshold) |
      (ifelse(is.na(dat$quakers),FALSE,dat$quakers) > js$defects$quakers$upper_threshold)
du[is.na(du)] <- FALSE
cat("union defect flag sum:", sum(du), " vs CSV:", sum(csv$defect_outlier_flag), "\n")
'

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")
A <- csv$altitude_outlier_flag; M <- csv$multivariate_grade_outlier_flag
D <- csv$defect_outlier_flag; I <- csv$high_leverage_flag & csv$high_cook_flag
tot <- A+M+D+I
fi <- js$flag_interactions
chk <- function(lab, a, b) cat(sprintf("  %-20s json=%-5d recompute=%-5d %s\n", lab, a, b, ifelse(a==b,"OK","MISMATCH")))
cat("Flag interactions vs CSV:\n")
chk("zero", fi$n_with_zero_flags, sum(tot==0)); chk("one", fi$n_with_one_flag, sum(tot==1))
chk("two", fi$n_with_two_flags, sum(tot==2)); chk("three", fi$n_with_three_flags, sum(tot==3))
chk("four", fi$n_with_four_flags, sum(tot==4))
chk("alt_only", fi$n_alt_only, sum(A&!M&!I&!D)); chk("mv_only", fi$n_mv_only, sum(M&!A&!I&!D))
chk("influence_only", fi$n_influence_only, sum(I&!A&!M&!D)); chk("defect_only", fi$n_defect_only, sum(D&!A&!M&!I))
chk("alt_and_mv", fi$n_alt_and_mv, sum(A&M)); chk("alt_and_influence", fi$n_alt_and_influence, sum(A&I))
chk("mv_and_influence", fi$n_mv_and_influence, sum(M&I)); chk("alt_mv_influence", fi$n_alt_mv_influence, sum(A&M&I))
cat("composite = A|M|(lev&cook):", all(csv$composite_drop_flag == (A|M|I)), " n_composite=", sum(csv$composite_drop_flag), "\n")
cat("influence n_both matches I:", js$influence$n_both == sum(I), "\n")
'

cd /app && cp outputs/outlier_report.json /tmp/run1.json && rm -rf outputs && mkdir -p outputs && Rscript analysis.R >/dev/null 2>&1 && diff <(cat /tmp/run1.json) <(cat outputs/outlier_report.json) && echo "REPRODUCIBLE: JSON identical across runs" || echo "DIFF FOUND"
echo "---sentinel check---"
Rscript -e '
suppressPackageStartupMessages(library(readr))
raw <- read_csv("/app/data/coffee_ratings.csv", show_col_types=FALSE)
g10 <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
az <- rowSums(raw[,g10]==0, na.rm=TRUE)==10
cat("all-grades-zero rows:", sum(az), " | their total_cup_points:", raw$total_cup_points[az], "\n")
cat("total_cup_points==0 rows:", sum(raw$total_cup_points==0), " | same row:", which(az)==which(raw$total_cup_points==0), "\n")
'

cd /app && Rscript -e '
fit <- readRDS("outputs/coffee_ols_selected.rds")
rn <- rownames(model.frame(fit))
cat("n model rows:", length(rn), "\n")
cat("first 5 rownames:", paste(head(rn,5), collapse=","), "\n")
cat("last 5 rownames:", paste(tail(rn,5), collapse=","), "\n")
cat("rownames == 1:1338 ?", all(rn == as.character(1:1338)), "\n")
cat("positional map: model row i -> row_id (i-1). Consistent since all 1338 complete & in order.\n")
# confirm response present for recompute
cat("has response total_cup_points in frame:", "total_cup_points" %in% names(model.frame(fit)), "\n")
cat("model stored (can recompute cooks.distance):", !is.null(fit$model), "\n")
'

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; ap <- alt[!is.na(alt)&alt>0]
q <- quantile(ap,c(.25,.75)); iqr<-q[[2]]-q[[1]]
cat("DRAFT raw-meter fence: [", round(q[[1]]-3*iqr,1), ",", round(q[[2]]+3*iqr,1), "] -> n flagged:", sum(!is.na(alt)&(alt<q[[1]]-3*iqr|alt>q[[2]]+3*iqr)),"\n")
cat("  (lower bound is negative/meaningless; high-altitude farms >3100m wrongly flagged, decimal slips like 1.0 caught but unfixable)\n")
# draft MD on 10 cols vs 7
g10<-c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
g7<-c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G10<-as.matrix(dat[,g10]); cc<-complete.cases(G10)
md10<-mahalanobis(G10[cc,],colMeans(G10[cc,]),cov(G10[cc,]),tol=1e-30)
cat("DRAFT MD(10 cols) flagged at chisq.999 df10:", sum(md10>qchisq(.999,10)), " vs 7-col:40\n")
# rstandard>2 count (draft influence)
fit<-lm(total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points,data=dat)
cat("DRAFT abs(rstandard)>2 count:", sum(abs(rstandard(fit))>2), " (~",round(100*sum(abs(rstandard(fit))>2)/nobs(fit),1),"% by construction, not influence)\n")
'

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

The draft applied one off-the-shelf recipe to every column. Each block below
says why that recipe breaks for the column's actual shape and what replaced it.
All figures come from the re-executed `analysis.R` (1,339 input lots).

## Data prep

One withdrawn submission has **all ten grades recorded as 0** (and
`total_cup_points = 0`). Left in, it is a leverage bomb for every model and a
fake extreme for every distance. I drop exactly that lot (identified by
all-grades-zero, equivalent here to the single `total_cup_points == 0` row):
**n_input = 1339 → n_after_zero_score_drop = 1338**. Survivors are numbered
`row_id` 0–1337 in input order so the flag table aligns 1:1 with the fits.

## Altitude

`altitude_mean_meters` is strongly right-skewed (median ≈ 1,311 m, max
190,164 m). The draft's Tukey fence on the **raw meter scale** gives
`[-400, 3100]`: the lower bound is physically impossible, genuine high-grown
lots above 3,100 m are flagged, and the decimal-slip errors are neither
isolated nor repaired. I build the k = 3 fence on `log10` of the positive
altitudes, then back-transform: **[357.4 m, 4923.8 m]**, flagging **51** lots.
For each flagged row I take the first numeric token of the raw `altitude`
string and test `÷10`, `÷100`, then as-is, keeping the first candidate inside
the fence. That recovers **7** unit slips (e.g. `190164 → 1901.64`,
`11000 metros → 1100`, `1100.00 mosl → 1100`); the rest (true lowland lots
near 1 m, range text like `1'500`) stay `NA`. The corrected meters feed the
altitude regression; the outlier flag is retained.

## Defects

`category_one_defects` (84.9% zero), `category_two_defects` (27.9% zero) and
`quakers` (93.0% zero) are mass-at-zero counts with a thin upper tail. The
draft's raw Tukey fence collapses: with Q1 = Q3 = 0 for the two sparse columns
the upper fence is **0**, so *every* nonzero count is "extreme" , **202** and
**94** lots (~15% and ~7%), which is meaningless. I instead cut each column at
its **98th percentile** and flag counts strictly above it: thresholds
5 / 20 / 2 flag **20 / 23 / 25** lots (1.5% / 1.7% / 1.9%), matching the
"top 1–2%" tail. A lot is a defect outlier if any column trips (**63** lots).

## Multivariate grades

The draft ran Mahalanobis on **all ten** grade columns. `uniformity`,
`clean_cup` and `sweetness` are near-constant at 10 (huge mass at the cap),
so they contribute almost no genuine covariance but inflate the covariance
condition number (≈ 67 vs ≈ 22 for the seven primaries) and manufacture
distance from cap deviations , 76 flags, many driven by clean-cup quirks. I
use a **classical** Mahalanobis distance (sample mean and covariance, *not*
MCD, as specified) on the **seven primary SCA attributes** only
(`aroma, flavor, aftertaste, acidity, body, balance, cupper_points`) with a
chi-squared 99.9% cutoff at df = 7 (**24.32**), flagging **40** lots.

## Influence

`abs(rstandard) > 2` measures large *residuals* in `y`, not influence on the
coefficients, and flags ~4% of rows (57) by construction regardless of
leverage. I first choose the model by AIC among three specs on the 1,338 rows
complete in the seven grades: grades_only **4459.4**, grades_altitude
**4461.3**, grades_defects **4418.7** → **grades_defects** wins (log1p defect
terms carry real signal; `total_cup_points` is the grade sum so the fit is a
near-identity plus defect adjustment). With n = 1338, p = 11 I flag **high
leverage** `h > 2p/n = 0.0164` (**110**) and **high Cook's D**
`> 4/n = 0.00299` (**66**) separately; the influence flag is the **conjunction**
(**35**), isolating points that are both unusual in X and actually move the
fit. The selected `lm` is saved so the verifier can recompute both diagnostics.

## Country ranking

The draft reported only raw means and never ran the robustness check it
promised (`trimmed_mean = raw_mean`), so every rank change was forced to 0. For
the **21** c
…[truncated]

cd /app && echo "Word count: $(wc -w < outputs/audit_memo.md)" && echo "---headings---" && grep '^## ' outputs/audit_memo.md && echo "---files in outputs---" && ls -la outputs/

cd /app && Rscript -e '
suppressPackageStartupMessages({library(jsonlite)})
js <- fromJSON("outputs/outlier_report.json")
pc <- js$country_ranking$per_country
pc$chg <- abs(pc$raw_rank - pc$trimmed_rank)
cat("Countries moving >=2 ranks:\n")
print(pc[pc$chg>=2, c("country","raw_rank","trimmed_rank","chg")], row.names=FALSE)
'

command: str_replace
path: /app/analysis.R
old_str: write_json(results, file.path(OUTPUT_DIR, "outlier_report.json"),
           auto_unbox = FALSE, pretty = TRUE, na = "null", digits = NA)
new_str: write_json(results, file.path(OUTPUT_DIR, "outlier_report.json"),
           auto_unbox = FALSE, pretty = TRUE, na = "null", digits = NA)

# ----------------------------------------------------------------------------
# 11. Audit memo (data-driven; regenerated on every run)
# ----------------------------------------------------------------------------
qd <- quantile(alt[!is.na(alt) & alt > 0], c(0.25, 0.75), names = FALSE)
draft_alt_lo <- qd[1] - 3 * (qd[2] - qd[1]); draft_alt_hi <- qd[2] + 3 * (qd[2] - qd[1])
draft_alt_n  <- sum(!is.na(alt) & (alt < draft_alt_lo | alt > draft_alt_hi))
draft_d1_n   <- sum(dat$category_one_defects > 0, na.rm = TRUE)
draft_q_n    <- sum(dat$quakers > 0, na.rm = TRUE)
G10 <- as.matrix(dat[, GRADES10]); cc10 <- complete.cases(G10)
md10 <- mahalanobis(G10[cc10, ], colMeans(G10[cc10, ]), cov(G10[cc10, ]), tol = 1e-30)
draft_md10_n <- sum(md10 > qchisq(0.999, df = length(GRADES10)))
draft_rstd_n <- sum(abs(rstandard(fit_grades_only)) > 2)
moving <- country_tbl %>% filter(abs_rank_change >= 2)
moving_txt <- paste(sprintf("%s (%d->%d)", moving$country, moving$raw_rank, moving$trimmed_rank),
                    collapse = ", ")
d1t <- d_list$category_one_defects$upper_threshold
d2t <- d_list$category_two_defects$upper_threshold
dqt <- d_list$quakers$upper_threshold
n_composite <- sum(composite)

memo <- c(
"# Coffee Quality Outlier Audit -- Method Memo",
"",
"The draft applied one off-the-shelf recipe to every column. Each block says",
"why that recipe breaks for the column's actual shape and what replaced it.",
"All figures are emitted by the re-executed `analysis.R`.",
"",
"## Data prep",
"",
sprintf("One withdrawn submission has **all ten grades recorded as 0** (and `total_cup_points = 0`). Left in, it is a leverage bomb for every model and a fake extreme for every distance. I drop exactly that lot (all-grades-zero, here the single `total_cup_points == 0` row): **n_input = %d -> n_after_zero_score_drop = %d**. Survivors are numbered `row_id` 0-%d in input order so the flag table aligns 1:1 with the fits.", n_input, n_after, n_after - 1L),
"",
"## Altitude",
"",
sprintf("`altitude_mean_meters` is strongly right-skewed (max 190,164 m). The draft's Tukey fence on the **raw meter scale** gives `[%.0f, %.0f]`: the lower bound is physically impossible, genuine high-grown lots are flagged, and decimal slips are never repaired (%d lots). I build the k = 3 fence on `log10` of the positive altitudes, then back-transform: **[%.1f m, %.1f m]**, flagging **%d** lots. For each flagged row I take the first numeric token of the raw `altitude` string and test `/10`, `/100`, then as-is, keeping the first candidate inside the fence. That recovers **%d** unit slips (e.g. 190164 -> 1901.64, '11000 metros' -> 1100, '1100.00 mosl' -> 1100); true lowland lots and range text stay `NA`. Corrected meters feed the altitude regression; the outlier flag is retained.", draft_alt_lo, draft_alt_hi, draft_alt_n, lo_m, hi_m, n_flag_alt, n_unit_corrected),
"",
"## Defects",
"",
sprintf("`category_one_defects` (84.9%% zero), `category_two_defects` (27.9%% zero) and `quakers` (93.0%% zero) are mass-at-zero counts with a thin upper tail. The draft's raw Tukey fence collapses: with Q1 = Q3 = 0 the upper fence is **0**, so *every* nonzero count is 'extreme' -- **%d** and **%d** lots, which is meaningless. I instead cut each column at its **98th percentile** and flag counts strictly above it: thresholds %g / %g / %g flag **%d / %d / %d** lots (~1.5-1.9%%), matching the top 1-2%% tail. A lot is a defect outlier if any column trips (**%d** lots).", draft_d1_n, draft_q_n, d1t, d2t, dqt, d_list$category_one_defects$n, d_list$category_two_defects$n, d_list$quakers$n, sum(flag_defect)),
"",
"## Multivariate grades",
"",
sprintf("The draft ran Mahalanob
…[truncated]

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R 2>&1 && echo "---FILES---" && ls outputs/ && echo "---MEMO WORDS---" && wc -w < outputs/audit_memo.md && echo "---HEADINGS---" && grep -c '^## ' outputs/audit_memo.md

cd /app && cat outputs/audit_memo.md

cd /app && cp outputs/outlier_report.json /tmp/r1.json && cp outputs/audit_memo.md /tmp/m1.md
rm -rf outputs && mkdir -p outputs && Rscript analysis.R >/dev/null 2>&1
echo "JSON reproducible: $(diff -q /tmp/r1.json outputs/outlier_report.json >/dev/null && echo YES || echo NO)"
echo "Memo reproducible: $(diff -q /tmp/m1.md outputs/audit_memo.md >/dev/null && echo YES || echo NO)"
echo "---final full reconciliation---"
Rscript -e '
suppressPackageStartupMessages({library(readr);library(jsonlite);library(dplyr)})
js<-fromJSON("outputs/outlier_report.json"); csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
fit<-readRDS("outputs/coffee_ols_selected.rds")
stopifnot(js$n_input==1339, js$n_after_zero_score_drop==1338, nrow(csv)==1338)
# schema key presence
need<-c("n_input","n_after_zero_score_drop","methods","altitude","defects","multivariate_grades","influence","country_ranking","sensitivity","flag_interactions")
cat("all top-level keys present:", all(need %in% names(js)),"\n")
cat("methods keys:", all(c("altitude","defects","multivariate_grades","influence","country_ranking") %in% names(js$methods)),"\n")
# model reconcile
n<-nobs(fit);p<-length(coef(fit))
lev<-hatvalues(fit)>2*p/n; ck<-cooks.distance(fit)>4/n
cat("lev/cook/composite reconcile:", all(lev==csv$high_leverage_flag)&&all(ck==csv$high_cook_flag)&&
    all(csv$composite_drop_flag==(csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag))),"\n")
# partition sum
fi<-js$flag_interactions
cat("flag partition sums to 1338:", fi$n_with_zero_flags+fi$n_with_one_flag+fi$n_with_two_flags+fi$n_with_three_flags+fi$n_with_four_flags==1338,"\n")
# per-country internal consistency: ranks match means
pc<-js$country_ranking$per_country
cat("raw_rank consistent with raw_mean order:", all(pc$raw_rank==rank(-pc$raw_mean,ties.method="min")),"\n")
cat("trimmed_rank consistent with trimmed_mean order:", all(pc$trimmed_rank==rank(-pc$trimmed_mean,ties.method="min")),"\n")
cat("n_countries==nrow(per_country):", js$country_ranking$n_countries_evaluated==nrow(pc),"\n")
# CI ordering and delta inside? (not required but sanity)
cat("CI lower<=upper:", js$sensitivity$delta_pct_ci_lower<=js$sensitivity$delta_pct_ci_upper,"\n")
cat("delta_pct recompute:", all.equal(js$sensitivity$delta_pct, 100*(js$sensitivity$coef_flavor_dropped-js$sensitivity$coef_flavor_full)/js$sensitivity$coef_flavor_full),"\n")
cat("selected_model valid:", js$influence$selected_model %in% c("grades_only","grades_altitude","grades_defects"),"\n")
cat("top10 length 10 & integer:", length(js$influence$top10_cook_row_ids)==10 && all(js$influence$top10_cook_row_ids==as.integer(js$influence$top10_cook_row_ids)),"\n")
cat("NA country rows in csv:", sum(is.na(csv$country_of_origin)),"\n")
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr);library(jsonlite);library(dplyr)})
js<-fromJSON("outputs/outlier_report.json")
# Re-derive bootstrap independently to confirm determinism & n=1338, reps=500
raw<-read_csv("/app/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),]
csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
comp<-csv$composite_drop_flag
f<-total_cup_points~aroma+flavor+aftertaste+acidity+body+balance+cupper_points
set.seed(20260512); B<-500; n<-nrow(dat); dl<-numeric(0)
for(b in 1:B){idx<-sample.int(n,n,replace=TRUE); s<-dat[idx,]; cb<-comp[idx]; k<-s[!cb,]
 if(nrow(k)<10) next
 ff<-tryCatch(lm(f,s),error=function(e)NULL); fd<-tryCatch(lm(f,k),error=function(e)NULL)
 if(is.null(ff)||is.null(fd)) next
 cf<-coef(ff)["flavor"]; cd<-coef(fd)["flavor"]; if(is.na(cf)||is.na(cd)||cf==0) next
 dl<-c(dl,100*(cd-cf)/cf)}
ci<-quantile(dl,c(.025,.975),names=FALSE)
cat("independent bootstrap: n=",n," reps=",length(dl)," CI=[",round(ci[1],4),",",round(ci[2],4),"]\n")
cat("matches JSON reps=500:",length(dl)==js$sensitivity$n_bootstrap_replicates,
    " CI match:",isTRUE(all.equal(ci[1],js$sensitivity$delta_pct_ci_lower))&&isTRUE(all.equal(ci[2],js$sensitivity$delta_pct_ci_upper)),"\n")
# defect semantics: confirm > (strict) matches reported counts and >= would differ
for(col in c("category_one_defects","category_two_defects","quakers")){
 t<-js$defects[[col]]$upper_threshold
 cat(sprintf("%-22s thr=%g  sum(x>thr)=%d (reported=%d)  sum(x>=thr)=%d\n",col,t,
   sum(dat[[col]]>t,na.rm=TRUE),js$defects[[col]]$n_flagged,sum(dat[[col]]>=t,na.rm=TRUE)))
}
'

cd /app && Rscript -e '
suppressPackageStartupMessages(library(jsonlite))
fit<-readRDS("outputs/coffee_ols_selected.rds"); js<-fromJSON("outputs/outlier_report.json")
ck<-cooks.distance(fit)
ids<-js$influence$top10_cook_row_ids
vals<-ck[ids+1]   # row_id = position-1
cat("top10 row_ids:     ", paste(ids,collapse=", "),"\n")
cat("their Cook D:      ", paste(round(vals,4),collapse=", "),"\n")
cat("strictly decreasing:", all(diff(vals)<0),"\n")
cat("these are the 10 largest overall:", setequal(ids, (order(ck,decreasing=TRUE)[1:10])-1),"\n")
'

cd /app && rm -rf outputs && mkdir -p outputs && Rscript -e 'options(warn=1); source("analysis.R")' 2>&1 | grep -iE "warn|error|cannot|fail" || echo "No warnings or errors."
echo "=== Deliverables in /app/outputs ==="
ls -la outputs/
echo "=== JSON validity ==="
Rscript -e '[redacted source](jsonlite::fromJSON("outputs/outlier_report.json")); cat("valid JSON\n")'

cd /app && echo "=== header + first 3 rows ===" && head -4 outputs/outlier_flags.csv
echo "=== a unit-corrected row ===" && awk -F, 'NR==1 || ($5!="" && $5!="NA")' outputs/outlier_flags.csv | head -3
echo "=== count of each flag (excl header) ===" 
Rscript -e '
suppressPackageStartupMessages(library(readr)); d<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
cat("rows:",nrow(d),"\n")
for(c in c("altitude_outlier_flag","defect_outlier_flag","multivariate_grade_outlier_flag","high_leverage_flag","high_cook_flag","composite_drop_flag")) cat(sprintf("  %-34s TRUE=%d\n",c,sum(d[[c]])))
cat("  altitude_corrected_m non-NA:",sum(!is.na(d$altitude_corrected_m)),"\n")
cat("  col types:",paste(sapply(d,function(x)class(x)[1]),collapse=", "),"\n")'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr);library(dplyr)})
d<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
raw<-read_csv("/app/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),]
corr<-d %>% filter(!is.na(altitude_corrected_m)) %>% select(row_id,altitude_corrected_m)
corr$raw_altitude <- dat$altitude[corr$row_id+1]
corr$orig_mean_m  <- dat$altitude_mean_meters[corr$row_id+1]
cat("Unit-corrected rows (all also altitude_outlier_flag=TRUE):\n")
print(as.data.frame(corr), row.names=FALSE)
cat("\nAll corrected rows flagged as altitude outliers:", all(d$altitude_outlier_flag[d$row_id %in% corr$row_id]),"\n")'

cd /app && Rscript -e '[redacted source](parse("analysis.R")); cat("analysis.R parses cleanly\n")'
echo "=== final end-to-end summary ==="
Rscript -e '
suppressPackageStartupMessages(library(jsonlite)); js<-fromJSON("outputs/outlier_report.json")
cat(sprintf("n_input=%d  n_after=%d\n",js$n_input,js$n_after_zero_score_drop))
cat(sprintf("altitude fence=[%.1f, %.1f]m  flagged=%d  corrected=%d\n",js$altitude$lower_bound_m,js$altitude$upper_bound_m,js$altitude$n_flagged,js$altitude$n_unit_corrected))
cat(sprintf("defects thresh/flag: c1=%g/%d c2=%g/%d q=%g/%d\n",js$defects$category_one_defects$upper_threshold,js$defects$category_one_defects$n_flagged,js$defects$category_two_defects$upper_threshold,js$defects$category_two_defects$n_flagged,js$defects$quakers$upper_threshold,js$defects$quakers$n_flagged))
cat(sprintf("MV: %d cols, thresh=%.2f, flagged=%d\n",length(js$multivariate_grades$columns_used),js$multivariate_grades$threshold,js$multivariate_grades$n_flagged))
cat(sprintf("influence: %s (AIC %.1f/%.1f/%.1f) lev=%d cook=%d both=%d\n",js$influence$selected_model,js$influence$aic_grades_only,js$influence$aic_grades_altitude,js$influence$aic_grades_defects,js$influence$n_high_leverage,js$influence$n_high_cook,js$influence$n_both))
cat(sprintf("country: %d eval, %d move>=2, max=%d\n",js$country_ranking$n_countries_evaluated,js$country_ranking$n_rank_changes_ge_2,js$country_ranking$max_abs_rank_change))
cat(sprintf("sensitivity: flavor %.3f->%.3f (%.1f%%) CI[%.1f,%.1f] reps=%d\n",js$sensitivity$coef_flavor_full,js$sensitivity$coef_flavor_dropped,js$sensitivity$delta_pct,js$sensitivity$delta_pct_ci_lower,js$sensitivity$delta_pct_ci_upper,js$sensitivity$n_bootstrap_replicates))
'