Pengukuran pengetahuan budidaya petani muda dengan Item Response Theory

# Contoh kode R: Item Response Theory (mirt)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(mirt)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis: 14 butir dikotomus (benar = 1, salah = 0) dari model 2PL ----
n <- 200 + 200                                                # IRT butuh sampel cukup besar
nb <- 14
z1 <- rnorm(n); z2 <- rnorm(n)
g1 <- runif(1, 0.25, 0.45); g2 <- runif(1, 0.1, 0.3) * sample(c(-1, 1), 1)
th <- g1 * z1 + g2 * z2 + rnorm(n, 0, sqrt(1 - g1^2 - g2^2))   # kemampuan laten (theta), varians 1
a_sb <- runif(nb, 0.7, 2.2)                                     # daya beda butir
b_sb <- sort(runif(nb, -2, 2))                                  # tingkat kesulitan butir
pr <- plogis(outer(th, b_sb, "-") * rep(a_sb, each = n))        # peluang menjawab benar
d <- data.frame(farming_experience = round(50 + 10 * z1, 2), years_of_schooling = round(30 + 6 * z2, 2))
jw <- matrix(rbinom(n * nb, 1, pr), n)
colnames(jw) <- sprintf("item%02d", 1:nb)
d <- cbind(d, jw)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Kalibrasi model 2PL ----
butir <- d[, sprintf("item%02d", 1:nb)]
fit <- mirt(butir, 1, itemtype = "2PL", verbose = FALSE)
par <- coef(fit, IRTpars = TRUE, simplify = TRUE)$items         # a = daya beda, b = kesulitan
print(round(par[, c("a", "b")], 3))

# ---- Kecocokan model dan butir ----
m2 <- M2(fit)                                                   # statistik M2 dan RMSEA
print(m2)
itf <- itemfit(fit, fit_stats = "S_X2")                         # S-X2 per butir
print(itf)
misfit <- sum(itf$p.S_X2 < 0.05, na.rm = TRUE)

# ---- Skor kemampuan (EAP) dan reliabilitas empiris ----
fs <- fscores(fit, method = "EAP", full.scores.SE = TRUE)
rel <- unname(empirical_rxx(fs))
d$cultivation_knowledge <- fs[, 1]                                              # theta hasil estimasi
reg <- lm(cultivation_knowledge ~ farming_experience + years_of_schooling, data = d)                    # kaitan theta dengan karakteristik peserta (b_x1 = per 1 SD X1)
print(summary(reg))
b_x1 <- coef(reg)[["farming_experience"]] * sd(d$farming_experience); p_x1 <- summary(reg)$coefficients["farming_experience", 4]


# ---- Grafik ----
tg <- seq(-4, 4, length.out = 121)
pb <- probtrace(fit, matrix(tg))                                # kolom: itemXX.P.0, itemXX.P.1
pk <- pb[, grep("P.1$", colnames(pb))]
dc <- data.frame(theta = rep(tg, nb), p = as.vector(pk), butir = rep(sprintf("item%02d", 1:nb), each = length(tg)),
                 kesulitan = rep(par[, "b"], each = length(tg)))
p1 <- ggplot(dc, aes(theta, p, group = butir, colour = kesulitan)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = warna$redup) +
  geom_line(linewidth = 0.75, alpha = 0.9) +
  scale_colour_gradient(low = pal[1], high = pal[2], name = "Difficulty b") +
  labs(title = "Item characteristic curves (2PL)", x = "cultivation_knowledge", y = "P(correct)")
ti <- seq(-3, 3, length.out = 91)
info <- testinfo(fit, matrix(ti))
di <- data.frame(theta = ti, info = info, se = 1 / sqrt(info))
k_se <- max(di$info) / max(di$se)                               # skala sumbu kedua untuk SE
p2 <- ggplot(di, aes(theta)) +
  geom_area(aes(y = info), fill = pal[3], alpha = 0.3) +
  geom_line(aes(y = info, colour = "Information"), linewidth = 1) +
  geom_line(aes(y = se * k_se, colour = "Standard error"), linewidth = 0.9, linetype = "dashed") +
  scale_colour_manual(values = c(Information = pal[3], `Standard error` = pal[4]), name = NULL) +
  scale_y_continuous(sec.axis = sec_axis(~ . / k_se, name = "SE")) +
  labs(title = "Test information and SE", x = "cultivation_knowledge", y = "Information")
p3 <- ggplot(d, aes(farming_experience, cultivation_knowledge)) +
  geom_point(colour = pal[1], alpha = 0.45, size = 1.4) +
  geom_smooth(method = "lm", formula = y ~ x, colour = pal[2], fill = pal[2], alpha = 0.18) +
  labs(title = "Estimated theta (EAP) vs covariate", x = "farming_experience", y = "Estimated theta")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3)

print(p_hasil)