Struktur dimensi kapasitas adaptasi petani terhadap perubahan iklim: EFA dan reliabilitas

# Contoh kode R: Exploratory Factor Analysis and Reliability
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(psych)
  library(GPArotation)
})
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: 12 butir Likert 1-5 yang dirancang untuk tiga dimensi ----
n <- 200
rho <- runif(3, 0.2, 0.45)                                      # korelasi antardimensi
R <- matrix(c(1, rho[1], rho[2], rho[1], 1, rho[3], rho[2], rho[3], 1), 3)
eta <- matrix(rnorm(n * 3), n) %*% chol(R)
likert <- function(e, k) {
  sapply(seq_len(k), function(j) {
    l <- runif(1, 0.62, 0.88)
    s <- l * e + rnorm(length(e), 0, sqrt(1 - l^2))
    findInterval(s, sort(c(-2, -1.1, -0.2, 0.9) + runif(1, -0.3, 0.2))) + 1
  })
}
d <- data.frame(likert(eta[, 1], 4), likert(eta[, 2], 4), likert(eta[, 3], 4))
names(d) <- c(paste0("climate_information_access", "_", 1:4), paste0("adaptation_resources", "_", 1:4), paste0("adaptation_skills", "_", 1:4))
d <- d[, sample(ncol(d))]                                       # urutan butir diacak seperti di kuesioner
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Kelayakan data: KMO dan uji Bartlett ----
km <- KMO(d)
bt <- cortest.bartlett(cor(d), n = nrow(d))
print(km); print(bt)

# ---- Jumlah faktor: analisis paralel ----
pa <- fa.parallel(d, fm = "ml", fa = "fa", plot = FALSE, n.iter = 20)
nf <- 3                                                         # sesuai rancangan instrumen (bandingkan dengan pa$nfact)

# ---- EFA: maximum likelihood, rotasi oblimin ----
ef <- fa(d, nfactors = nf, rotate = "oblimin", fm = "ml")
print(ef$loadings, cutoff = 0.3, sort = TRUE)
L <- unclass(ef$loadings)
dim_butir <- sub("_[0-9]+$", "", names(d))                      # dimensi rancangan tiap butir
dims <- c("climate_information_access", "adaptation_resources", "adaptation_skills")
# cocokkan faktor hasil EFA dengan dimensi (muatan rata-rata terbesar)
urut <- sapply(dims, function(g) which.max(colMeans(abs(L[dim_butir == g, , drop = FALSE]))))
if (length(unique(urut)) < 3) urut <- 1:3
L <- L[, urut]; colnames(L) <- dims
benar <- mean(apply(abs(L), 1, which.max) == match(dim_butir, dims))   # proporsi butir di faktor yang benar
varjel <- sum(colSums(L^2)) / ncol(d)                           # proporsi varians dijelaskan tiga faktor

# ---- Reliabilitas: alpha Cronbach per dimensi ----
alfa <- sapply(dims, function(g) psych::alpha(d[, dim_butir == g], warnings = FALSE)$total$raw_alpha)
print(round(alfa, 3))


# ---- Grafik ----
singkat <- function(s, w = 18) {                                # nama dimensi pendek untuk grafik
  k <- strsplit(s, " ")[[1]]
  if (nchar(s) > w) k <- k[!tolower(k) %in% c("of", "to", "the", "and", "for", "in", "on", "with")]
  i <- length(k)
  while (nchar(paste(k, collapse = " ")) > w && i >= 1) {
    if (nchar(k[i]) > 5) k[i] <- paste0(substr(k[i], 1, 4), ".")
    i <- i - 1
  }
  substr(paste(k, collapse = " "), 1, w)
}
nm_k <- make.unique(unname(sapply(c("climate_information_access", "adaptation_resources", "adaptation_skills"), singkat)), sep = " ")
ds <- data.frame(faktor = rep(seq_along(pa$fa.values), 2), nilai = c(pa$fa.values, pa$fa.sim),
                 jenis = rep(c("Observed", "Simulated (parallel)"), each = length(pa$fa.values)))
p1 <- ggplot(ds, aes(faktor, nilai, colour = jenis, linetype = jenis)) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_line(linewidth = 0.9) + geom_point(size = 2) +
  annotate("text", x = pa$nfact + 0.3, y = pa$fa.values[pa$nfact], label = paste(pa$nfact, "factors"),
           hjust = 0, colour = warna$teks, size = 3.2) +
  scale_colour_manual(values = pal[c(1, 2)], name = NULL) +
  scale_linetype_manual(values = c("solid", "dashed"), name = NULL) +
  scale_x_continuous(breaks = seq_along(pa$fa.values)) +
  labs(title = "Scree plot with parallel analysis", x = "Factor number", y = "Eigenvalue (FA)")
kode <- paste0("D", match(dim_butir, dims), ".", sub(".*_", "", names(d)))   # kode butir, mis. D2.3
dl <- data.frame(butir = rep(kode, 3), faktor = rep(paste0("F", 1:3), each = ncol(d)),
                 muatan = as.vector(L))
dl$butir <- factor(dl$butir, levels = rev(sort(unique(kode))))
p2 <- ggplot(dl, aes(faktor, butir, fill = muatan)) +
  geom_tile(colour = warna$panel, linewidth = 0.8) +
  geom_text(aes(label = sprintf("%.2f", muatan), fontface = ifelse(abs(muatan) >= 0.4, "bold", "plain")),
            colour = ifelse(abs(dl$muatan) >= 0.4, warna$panel, warna$teks), size = 2.8) +
  scale_fill_gradient2(low = pal[2], mid = warna$panel, high = pal[1], midpoint = 0, limits = c(-1, 1), guide = "none") +
  labs(title = "Pattern loadings (oblimin)", x = NULL, y = NULL,
       subtitle = paste(strwrap(paste0("D", 1:3, "/F", 1:3, " = ", nm_k, collapse = "; "), 48), collapse = "\n"))
da <- data.frame(dimensi = factor(nm_k, levels = nm_k), alfa = alfa)
p3 <- ggplot(da, aes(dimensi, alfa, fill = dimensi)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  geom_hline(yintercept = 0.7, linetype = "dashed", colour = warna$teks) +
  geom_text(aes(label = sprintf("%.3f", alfa)), vjust = -0.4, size = 3.1, colour = warna$teks) +
  scale_fill_manual(values = pal[1:3]) +
  scale_y_continuous(limits = c(0, 1.05), breaks = seq(0, 1, 0.25)) +
  labs(title = "Cronbach alpha per dimension", x = NULL, y = "Alpha")
p_hasil <- if (tata == 1) (p1 / p3) | p2 else (p1 | p2) / p3 + plot_layout(heights = c(1.5, 1))

print(p_hasil)