1 Sumber data

Analisis ini menggunakan data kualitas kopi yang bersumber dari Coffee Quality Institute (CQI), lembaga sertifikasi mutu kopi internasional. Setiap baris data merepresentasikan satu sampel kopi yang dinilai melalui proses cupping oleh panel Q-grader bersertifikat. Variabel yang tercatat meliputi skor sensorik (aroma, rasa, keasaman, body, dan lainnya), karakteristik asal (negara produsen, ketinggian lahan, metode pengolahan), serta jumlah cacat biji. Data dikompilasi oleh James LeDoux dari basis data resmi CQI dan dipublikasikan secara terbuka melalui proyek TidyTuesday, sehingga tersedia dalam format CSV yang dapat diunduh dan diolah langsung di R.

URL <- "https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-07-07/coffee_ratings.csv"
d <- read.csv(URL, stringsAsFactors = FALSE)
cat("Jumlah baris mentah:", nrow(d))
## Jumlah baris mentah: 1339

1.1 Pembersihan & penyiapan variabel

d <- subset(d, total_cup_points > 0)                       # buang skor 0 (galat input)
d$altitude_mean_meters[d$altitude_mean_meters < 1 |
                       d$altitude_mean_meters > 3000] <- NA # altitude tak masuk akal -> NA
d$alt100     <- d$altitude_mean_meters / 100               # altitude per 100 meter
d$moist_pct  <- d$moisture * 100                           # kelembapan dalam persen

# Negara: 5 produsen terbanyak + "Other" (referensi = Mexico)
top <- names(sort(table(d$country_of_origin), decreasing = TRUE))[1:5]
d$country <- factor(ifelse(d$country_of_origin %in% top, d$country_of_origin, "Other"))
d$country <- relevel(d$country, ref = "Mexico")

# Metode proses (respons NOMINAL): 3 level utama
keep <- c("Washed / Wet", "Natural / Dry", "Semi-washed / Semi-pulped")
d$proc <- factor(ifelse(d$processing_method %in% keep, d$processing_method, NA),
                 levels = keep)

# Grade kualitas (respons ORDINAL) -> ambang mengacu skala mutu kopi (skor 0-100)
d$grade <- factor(cut(d$total_cup_points, c(-Inf, 82, 84, Inf),
              labels = c("Standard (<82)", "Premium (82-83.9)", "Excellent (>=84)")),
              ordered = TRUE)
Model Respons Tipe Variabel
Multinomial metode proses Nominal (3) Washed/Wet, Natural/Dry, Semi-washed
Ordinal grade mutu Ordinal (3) Standard < Premium < Excellent
Poisson jumlah cacat Count category_two_defects

2 Regresi Logistik Multinomial — metode pengolahan biji

Pertanyaan: apakah ketinggian lahan, kelembapan, dan jumlah cacat memengaruhi metode pengolahan yang dipakai? Respons nominal; referensi = Washed / Wet.

Catatan pemodelan: country sengaja tidak dipakai di model multinomial karena negara hampir menentukan metode proses (mis. Brasil ⇒ natural), sehingga timbul complete separation yang membuat galat baku meledak. Prediktor kontinu dipakai agar estimasi stabil.

a <- droplevels(na.omit(d[, c("proc", "alt100", "moist_pct", "category_two_defects")]))
a$proc <- relevel(a$proc, ref = "Washed / Wet")

m1 <- multinom(proc ~ alt100 + moist_pct + category_two_defects, data = a, trace = FALSE)
co <- summary(m1)$coefficients; se <- summary(m1)$standard.errors
p  <- 2 * (1 - pnorm(abs(co / se)))

tab1 <- data.frame(
  Kontras  = rep(rownames(co), each = ncol(co)),
  Variabel = rep(colnames(co), times = nrow(co)),
  b        = round(as.vector(t(co)), 4),
  RRR      = round(exp(as.vector(t(co))), 3),
  p_value  = round(as.vector(t(p)), 4))
kable(tab1, caption = sprintf("Multinomial (n = %d, ref outcome = Washed/Wet)", nrow(a)))
Multinomial (n = 963, ref outcome = Washed/Wet)
Kontras Variabel b RRR p_value
Natural / Dry (Intercept) 0.7272 2.069 0.0146
Natural / Dry alt100 -0.1185 0.888 0.0000
Natural / Dry moist_pct -0.0705 0.932 0.0001
Natural / Dry category_two_defects 0.0057 1.006 0.7116
Semi-washed / Semi-pulped (Intercept) -1.3786 0.252 0.0161
Semi-washed / Semi-pulped alt100 -0.1456 0.864 0.0000
Semi-washed / Semi-pulped moist_pct 0.0609 1.063 0.1564
Semi-washed / Semi-pulped category_two_defects -0.0213 0.979 0.4632

Interpretasi. Tiap kenaikan ketinggian 100 m menurunkan kecenderungan biji diproses Natural/Dry (RRR 0,89; p < 0,001) maupun Semi-washed (RRR 0,86; p < 0,001) dibanding Washed/Wet — artinya kopi dataran tinggi condong diproses basah (washed), sedangkan natural/dry lebih lazim di lahan lebih rendah dan berkelembapan lebih rendah (RRR 0,93 per 1% kelembapan; p < 0,001). Jumlah cacat tidak berkaitan signifikan dengan pilihan proses.

3 Regresi Logistik Ordinal — grade mutu kopi

Pertanyaan: apa yang mendorong kopi naik ke grade mutu lebih tinggi? Respons berurut: Standard < Premium < Excellent. Koefisien positif ⇒ peluang lebih besar ke grade lebih tinggi.

b <- droplevels(na.omit(d[, c("grade", "proc", "alt100", "country", "moist_pct")]))
m2 <- polr(grade ~ proc + alt100 + country + moist_pct, data = b, Hess = TRUE)
ct <- coef(summary(m2))
ct <- cbind(ct, p = 2 * (1 - pnorm(abs(ct[, "t value"]))), OR = exp(ct[, "Value"]))
kable(round(ct, 4),
      caption = sprintf("Ordinal proportional-odds (n = %d). Dua baris akhir = ambang.", nrow(b)))
Ordinal proportional-odds (n = 963). Dua baris akhir = ambang.
Value Std. Error t value p OR
procNatural / Dry 0.4970 0.1867 2.6623 0.0078 1.6438
procSemi-washed / Semi-pulped 0.4533 0.3100 1.4623 0.1437 1.5735
alt100 0.1094 0.0188 5.8216 0.0000 1.1156
countryBrazil 0.8620 0.2838 3.0371 0.0024 2.3680
countryColombia 1.4983 0.2356 6.3601 0.0000 4.4740
countryGuatemala 0.6796 0.2120 3.2059 0.0013 1.9731
countryOther 0.9745 0.1823 5.3446 0.0000 2.6499
countryTaiwan 0.6432 0.3204 2.0079 0.0447 1.9026
moist_pct -0.0207 0.0155 -1.3428 0.1794 0.9795
Standard (<82)|Premium (82-83.9) 1.6282 0.3201 5.0861 0.0000 5.0945
Premium (82-83.9)|Excellent (>=84) 3.8829 0.3424 11.3416 0.0000 48.5639

Interpretasi. Ketinggian adalah pendorong mutu yang kuat: tiap +100 m menaikkan odds naik grade ~12% (OR 1,12; p < 0,001) — sesuai pengetahuan umum bahwa kopi dataran tinggi cenderung lebih bermutu. Pengolahan Natural/Dry berasosiasi dengan grade lebih tinggi (OR 1,64; p = 0,008). Asal negara sangat berpengaruh dibanding Meksiko (referensi): Kolombia OR 4,47; “Other” 2,65; Brasil 2,37; Guatemala 1,97; Taiwan 1,90 (semua signifikan). Kelembapan tidak signifikan.

4 Regresi Poisson — jumlah cacat biji

Pertanyaan: faktor apa yang menaikkan jumlah cacat tipe-2 (category_two_defects, data cacah)?

cc <- droplevels(na.omit(d[, c("category_two_defects", "proc", "alt100",
                               "country", "moist_pct")]))
c(mean = mean(cc$category_two_defects), varians = var(cc$category_two_defects))
##      mean   varians 
##  3.719626 29.734199
m3 <- glm(category_two_defects ~ proc + alt100 + country + moist_pct,
          data = cc, family = poisson)
s3 <- summary(m3)$coefficients
kable(round(cbind(s3, IRR = exp(s3[, "Estimate"])), 4),
      caption = sprintf("Poisson (n = %d). IRR = Incidence Rate Ratio.", nrow(cc)))
Poisson (n = 963). IRR = Incidence Rate Ratio.
Estimate Std. Error z value Pr(>|z|) IRR
(Intercept) 1.8399 0.0876 21.0028 0.0000 6.2956
procNatural / Dry 0.2636 0.0491 5.3703 0.0000 1.3016
procSemi-washed / Semi-pulped 0.0054 0.0818 0.0664 0.9471 1.0054
alt100 -0.0101 0.0050 -2.0338 0.0420 0.9899
countryBrazil -1.0456 0.0795 -13.1595 0.0000 0.3515
countryColombia -1.2714 0.0772 -16.4666 0.0000 0.2804
countryGuatemala -0.7194 0.0542 -13.2724 0.0000 0.4870
countryOther -0.7236 0.0427 -16.9658 0.0000 0.4850
countryTaiwan -2.8487 0.2084 -13.6665 0.0000 0.0579
moist_pct 0.0170 0.0049 3.4405 0.0006 1.0171 
disp <- deviance(m3) / df.residual(m3)
m3nb <- glm.nb(category_two_defects ~ proc + alt100 + country + moist_pct, data = cc)
data.frame(dispersi_dev_df = round(disp, 2),
           AIC_Poisson     = round(AIC(m3), 1),
           AIC_NegBin      = round(AIC(m3nb), 1))
dispersi_dev_df AIC_Poisson AIC_NegBin
4.52 6526.3 4475.4

Interpretasi. Pengolahan Natural/Dry menaikkan laju cacat ~30% (IRR 1,30; p < 0,001) dibanding washed — wajar karena pengeringan dengan kulit buah lebih rawan cacat. Tiap +100 m ketinggian sedikit menurunkan cacat (IRR 0,99; p = 0,04), dan tiap +1% kelembapan menaikkan cacat ~1,7% (p < 0,001). Dibanding Meksiko, semua negara lain punya laju cacat jauh lebih rendah (mis. Taiwan IRR 0,06; Kolombia 0,28).

Diagnostik overdispersi. Varians cacat (≈30) jauh di atas rata-ratanya (≈3,7), dan dispersi deviance/df ≈ 4.5 (idealnya ≈1). Asumsi Poisson dilanggar, sehingga galat baku Poisson terlalu optimis. Negative Binomial jauh lebih cocok (AIC 4475 vs Poisson 6526). Untuk data cacah seperti ini, gunakan Negative Binomial / quasi-Poisson.

5 Ringkasan

  • Multinomial (nnet::multinom) untuk respons nominal — metode proses ditentukan terutama oleh ketinggian & kelembapan.
  • Ordinal (MASS::polr) untuk respons berurut — mutu kopi naik seiring ketinggian, pengolahan natural, dan asal negara tertentu.
  • Poisson (glm, lalu MASS::glm.nb) untuk count — cacat dipengaruhi proses/asal/kelembapan; data overdispersi ⇒ pakai Negative Binomial.

6 Cara publikasi ke RPubs

  1. Buka file .Rmd ini di RStudio (paket: rmarkdown, knitr, nnet, MASS).
  2. Klik Knit → muncul HTML.
  3. Di viewer, klik Publish → pilih RPubs → beri judul → Publish. RPubs akan memberi URL publik gratis untuk dibagikan.