Peringkasan indikator mutu sekolah dengan PCA

# Contoh kode R: Principal Component Analysis (PCA)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(psych)
  library(ggrepel)
})
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: 7 sampai 9 indikator berkorelasi dari dua dimensi laten ----
n <- 200
positif <- function(x, m) if (m > 0) pmax(x, 0.05 * m) else x   # cegah nilai negatif yang tidak wajar
p_tambah <- sample(4:6, 1)                       # jumlah indikator tambahan
f1 <- rnorm(n); f2 <- 0.25 * f1 + rnorm(n)        # dua dimensi laten yang sedikit berkorelasi
indikator <- function(l1, l2) {
  # skor z indikator: muatan pada f1 dan f2 ditambah galat unik
  u <- sqrt(pmax(1 - l1^2 - l2^2, 0.15))
  l1 * f1 + l2 * f2 + rnorm(n, 0, u)
}
tanda <- sample(c(-1, 1), 3, replace = TRUE, prob = c(0.3, 0.7))
d <- data.frame(accreditation_score = positif(100 + 15 * indikator(tanda[1] * runif(1, 0.65, 0.85), runif(1, 0, 0.2)), 100),
                degree_qualified_teachers = positif(50 + 10 * indikator(tanda[2] * runif(1, 0.6, 0.8), runif(1, 0, 0.25)), 50),
                learning_facilities_index = positif(30 + 6 * indikator(runif(1, 0, 0.25), tanda[3] * runif(1, 0.6, 0.8)), 30))
for (j in seq_len(p_tambah)) {
  ke_f1 <- j %% 2 == 1
  l1 <- if (ke_f1) runif(1, 0.55, 0.85) * sample(c(-1, 1), 1, prob = c(0.25, 0.75)) else runif(1, 0, 0.25)
  l2 <- if (ke_f1) runif(1, 0, 0.25) else runif(1, 0.55, 0.85) * sample(c(-1, 1), 1, prob = c(0.25, 0.75))
  d[[sprintf("ind%02d", j + 3)]] <- round(50 + 10 * indikator(l1, l2), 2)
}
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Kelayakan data: KMO dan uji Bartlett ----
r <- cor(d)
kmo <- psych::KMO(r)
bart <- psych::cortest.bartlett(r, n = nrow(d))
print(kmo)
print(bart)

# ---- PCA pada data terstandar (matriks korelasi) ----
pca <- prcomp(d, center = TRUE, scale. = TRUE)
for (j in seq_len(ncol(pca$rotation))) if (sum(pca$rotation[, j]) < 0) {   # tanda komponen dibuat searah mayoritas indikator
  pca$rotation[, j] <- -pca$rotation[, j]; pca$x[, j] <- -pca$x[, j]
}
print(summary(pca))
eig <- pca$sdev^2
var_pct <- 100 * eig / sum(eig)
kum_pct <- cumsum(var_pct)
m_kaiser <- max(1, sum(eig > 1))                 # kriteria Kaiser (eigenvalue > 1)
muatan <- pca$rotation %*% diag(pca$sdev)        # loading = korelasi variabel dengan komponen
colnames(muatan) <- paste0("PC", seq_len(ncol(muatan)))
print(round(muatan[, seq_len(max(2, m_kaiser))], 3))
top_pc1 <- rownames(muatan)[which.max(abs(muatan[, 1]))]


# ---- Grafik ----
label_var <- function(v) {
  # nama variabel pendek untuk grafik (dibungkus per 14 karakter)
  nm <- c(accreditation_score = "accreditation_score", degree_qualified_teachers = "degree_qualified_teachers", learning_facilities_index = "learning_facilities_index")
  s <- ifelse(v %in% names(nm), nm[v], v)
  sapply(s, function(x) paste(strwrap(x, 14), collapse = "\n"))
}
ds <- data.frame(pc = factor(seq_along(eig)), var = var_pct, kum = kum_pct, simpan = seq_along(eig) <= m_kaiser)
p1 <- ggplot(ds, aes(pc, var)) +
  geom_col(aes(fill = simpan), width = 0.7, alpha = 0.85) +
  geom_line(aes(y = kum, group = 1), colour = pal[2], linewidth = 0.9) +
  geom_point(aes(y = kum), colour = pal[2], size = 1.8) +
  geom_hline(yintercept = 100 / length(eig), linetype = "dashed", colour = warna$redup) +
  scale_fill_manual(values = c(`TRUE` = pal[1], `FALSE` = warna$redup), guide = "none") +
  scale_y_continuous(limits = c(0, 105)) +
  labs(title = "Scree plot and cumulative variance", subtitle = "Bars = % variance, line = cumulative %, dashed = eigenvalue 1",
       x = "Component", y = "Variance explained (%)")
sk <- as.data.frame(pca$x[, 1:2])
ld <- data.frame(var = rownames(muatan), PC1 = muatan[, 1], PC2 = muatan[, 2])
ld$label <- label_var(ld$var)
f_skala <- 0.85 * max(abs(sk)) / max(abs(ld[, c("PC1", "PC2")]))
ld$x <- ld$PC1 * f_skala; ld$y <- ld$PC2 * f_skala
p2 <- ggplot() +
  geom_hline(yintercept = 0, colour = warna$grid) + geom_vline(xintercept = 0, colour = warna$grid) +
  geom_point(data = sk, aes(PC1, PC2), colour = pal[3], alpha = 0.4, size = 1.3) +
  geom_segment(data = ld, aes(x = 0, y = 0, xend = x, yend = y), colour = pal[2], linewidth = 0.8,
               arrow = arrow(length = unit(0.18, "cm"))) +
  ggrepel::geom_text_repel(data = ld, aes(x, y, label = label), colour = warna$teks, size = 2.8, lineheight = 0.85,
                           seed = 1, box.padding = 0.35, min.segment.length = 0.3, segment.colour = warna$redup) +
  labs(title = "PCA biplot (scores and loadings)", x = sprintf("PC1 (%.1f%%)", var_pct[1]), y = sprintf("PC2 (%.1f%%)", var_pct[2]))
m_tampil <- min(3, max(2, m_kaiser), ncol(muatan))
lm_ <- data.frame(var = rep(rownames(muatan), m_tampil), pc = rep(colnames(muatan)[seq_len(m_tampil)], each = nrow(muatan)),
                  nilai = as.vector(muatan[, seq_len(m_tampil)]))
lm_$label <- factor(label_var(lm_$var), levels = rev(unique(label_var(rownames(muatan)))))
p3 <- ggplot(lm_, aes(pc, label, fill = nilai)) +
  geom_tile(colour = warna$panel, linewidth = 0.8) +
  geom_text(aes(label = sprintf("%.2f", round(nilai, 2) + 0)), colour = warna$teks, size = 2.8) +
  scale_fill_gradient2(low = pal[1], mid = warna$panel, high = pal[2], midpoint = 0, limits = c(-1, 1), guide = "none") +
  labs(title = "Loadings (correlation with PC)", x = NULL, y = NULL) +
  theme(panel.grid = element_blank(), axis.text.y = element_text(size = rel(0.82), lineheight = 0.85))
p_hasil <- if (tata == 1) (p2 | (p1 / p3)) + plot_layout(widths = c(1.4, 1)) else
  ((p1 | p3) / p2) + plot_layout(heights = c(1, 1.3))

print(p_hasil)