Import Data & Seleksi Variabel

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(psych)
## 
## Attaching package: 'psych'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(corrplot)
## corrplot 0.95 loaded
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:psych':
## 
##     logit
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(FactoMineR)
library(GGally)
library(dplyr)
library(readr)
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(glmnet)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loaded glmnet 4.1-10
library(Metrics)
## 
## Attaching package: 'Metrics'
## 
## The following objects are masked from 'package:caret':
## 
##     precision, recall
library(ggcorrplot)
data <- read.csv("Instagram_Analytics.csv")
vars <- c("follower_count","likes","comments","shares","saves",
            "reach","impressions","engagement_rate","followers_gained","hashtags_count")
df <- data[vars]

Cek Missing Value

colSums(is.na(df))
##   follower_count            likes         comments           shares 
##                0                0                0                0 
##            saves            reach      impressions  engagement_rate 
##                0                0                0                0 
## followers_gained   hashtags_count 
##                0                0

Statistik Deskriptif

df %>% 
                 gather(variable, value) %>% 
                 ggplot(aes(x=value)) +
                 facet_wrap(~variable, scales="free") +
                 geom_histogram(bins=30, fill="steelblue", color="white")

Visualisasi histogram menunjukkan bahwa sebagian besar variabel, terutama likes, comments, shares, saves, reach, dan impressions, memiliki distribusi yang condong ke kanan (positively skewed), menandakan adanya sejumlah kecil postingan dengan performa sangat tinggi dibandingkan mayoritas lainnya. Sementara itu, engagement_rate dan hashtags_count relatif lebih terkonsentrasi, dan followers_gained cenderung lebih merata. Perbedaan bentuk distribusi serta skala antarvariabel ini mengindikasikan perlunya proses standardisasi sebelum analisis PCA agar tidak terjadi dominasi variabel dengan variansi besar.

Uji Asumsi dan Kelayakan Analisis Faktor

Rename Variabel (X1-X10)

orig_names <- colnames(df)
colnames(df) <- paste0("X", seq_len(ncol(df)))
mapping <- data.frame(original = orig_names, new = colnames(df))
print(mapping)
##            original new
## 1    follower_count  X1
## 2             likes  X2
## 3          comments  X3
## 4            shares  X4
## 5             saves  X5
## 6             reach  X6
## 7       impressions  X7
## 8   engagement_rate  X8
## 9  followers_gained  X9
## 10   hashtags_count X10

Matriks Korelasi

cor_mat_X <- cor(df, use = "pairwise.complete.obs", method = "pearson")
round(cor_mat_X, 3)
##         X1     X2     X3     X4     X5     X6     X7     X8     X9    X10
## X1   1.000 -0.002 -0.003 -0.005 -0.003  0.002  0.001 -0.006  0.013  0.006
## X2  -0.002  1.000  0.910  0.945  0.948  0.738  0.723  0.506 -0.005  0.006
## X3  -0.003  0.910  1.000  0.886  0.886  0.690  0.674  0.452 -0.006  0.004
## X4  -0.005  0.945  0.886  1.000  0.921  0.716  0.700  0.469 -0.003  0.007
## X5  -0.003  0.948  0.886  0.921  1.000  0.722  0.707  0.478 -0.003  0.009
## X6   0.002  0.738  0.690  0.716  0.722  1.000  0.986  0.006 -0.005  0.002
## X7   0.001  0.723  0.674  0.700  0.707  0.986  1.000 -0.021 -0.006  0.002
## X8  -0.006  0.506  0.452  0.469  0.478  0.006 -0.021  1.000  0.001 -0.001
## X9   0.013 -0.005 -0.006 -0.003 -0.003 -0.005 -0.006  0.001  1.000 -0.004
## X10  0.006  0.006  0.004  0.007  0.009  0.002  0.002 -0.001 -0.004  1.000
if (exists("cor_mat_X")) {
       cor_mat <- cor_mat_X
   } else {
           cor_mat <- cor(df, use = "pairwise.complete.obs", method = "pearson")
       }
if (exists("orig_names")) {
       colnames(cor_mat) <- rownames(cor_mat) <- orig_names
   }
corrplot(cor_mat, 
                     method = "color", 
                     type = "full", 
                     tl.cex = 0.8, 
                    tl.col = "black", 
                     addCoef.col = "black", 
                     col = colorRampPalette(c("blue", "white", "red"))(200))

Berdasarkan matriks korelasi ini r = 0,986 menunjukkan hampir hubungan linear sempurna dengan reach dan impressions yang hampir identik. Kemudian, Likes sangat kuat berkorelasi dengan saves, shares, dan comments yang membentuk satu klaster engagement. Hal ini menunjukkan indikasi yang kuat akan terbentuk faktor engagement dan faktor exposure dalam PCA/FA.

Uji Multikolinearitas (VIF Awal)

df_num <- df %>% select(where(is.numeric))
vif_vals <- sapply(names(df_num), function(v){
       R2 <- summary(lm(as.formula(paste(v, "~ .")), data = df_num))$r.squared
       ifelse(is.na(R2), NA, 1 / (1 - R2))
   })
vif_table <- data.frame(variable = names(vif_vals), VIF = round(vif_vals, 3))
print(vif_table)
##     variable    VIF
## X1        X1  1.000
## X2        X2 19.021
## X3        X3  6.171
## X4        X4 10.174
## X5        X5 10.806
## X6        X6 36.742
## X7        X7 36.151
## X8        X8  2.334
## X9        X9  1.000
## X10      X10  1.000

Nilai vif yang lebih dari 10 menandakan bahwa multikolinearitas sangat tinggi. Berdasarkan tabel vif di atas, nilai Reach & impressions sangat ekstrem karena r = 0,986. Begitu juga dengan Likes memiliki nilai vif yang tinggi karena korelasi kuat dengan 3 variabel lain. Sehingga, hanya nilai VIF = 1 terdapat pada X1, X9, X10 yang tidak berkorelasi dengan variabel lain. Dengan demikian, strukturnya sangat redundan, sehingga perlu penanganan yang lebih lanjut.

Log Transformation pada VIF Awal

df_log <- df %>% mutate(across(everything(), ~log1p(.)))
df_vif <- df_log %>% select(where(is.numeric))
vif_vals <- sapply(names(df_vif), function(v){
       R2 <- summary(lm(as.formula(paste(v, "~ .")), data = df_vif))$r.squared
       ifelse(is.na(R2), NA, 1/(1-R2))
   })
print(round(vif_vals,3))
##     X1     X2     X3     X4     X5     X6     X7     X8     X9    X10 
##  1.000  3.129  2.680  3.605  2.519 44.746 49.849  4.206  1.001  1.000

Setelah transformasi log, sebagian besar variabel menunjukkan penurunan multikolinearitas dengan nilai vif yang berada pada kisaran 2–4 sehingga tergolong aman. Namun, pada X6 dan X7 tetap memiliki vif yang sangat tinggi (44,746 dan 49,849), yang menunjukkan adanya hubungan linear hampir sempurna dan redundansi informasi. Sementara itu, X1, X9, dan X10 memiliki nilai vif ≈ 1 yang menandakan tidak terdapat korelasi linear yang berarti. Dengan demikian, transformasi log hanya efektif mengurangi multikolinearitas pada sebagian variabel, tetapi tidak pada X6 dan X7.

Residualisasi X7 terhadap X6

fit_log <- lm(X7 ~ X6, data = df_log)
df_log$X7_resid <- resid(fit_log)
df_vif10_log <- df_log %>% select(X1, X2, X3, X4, X5, X6, X7_resid, X8, X9, X10) #
dummy <- rnorm(nrow(df_vif10_log))
vif_vals_10_log <- vif(lm(dummy ~ ., data = df_vif10_log))
round(vif_vals_10_log, 3)
##       X1       X2       X3       X4       X5       X6 X7_resid       X8 
##    1.000    3.129    2.680    3.605    2.519    5.613    1.158    4.206 
##       X9      X10 
##    1.001    1.000

Setelah dilakukan residualisasi X7 terhadap X6 pada ruang log, nilai vif menunjukkan perbaikan signifikan pada masalah multikolinearitas. Variabel X7_resid memiliki vif sebesar 1,158 yang menandakan tidak lagi yang berkorelasi kuat dengan X6, sementara nilai vif X6 turun menjadi 5,613 (mendekati batas moderat). Variabel lainnya tetap berada pada kisaran aman (sekitar 1–4). Hal ini menunjukkan bahwa teknik residualisasi berhasil menghilangkan redundansi informasi antara X6 dan X7, sehingga struktur data menjadi lebih stabil dan layak untuk analisis lanjutan pada PCA dan analisis faktor.

Uji KMO

vars_kmo <- c("X1","X2","X3","X4","X5","X6","X7_resid","X8","X9","X10")
mat <- cor(df_log[vars_kmo], use = "pairwise.complete.obs")
kmo_res <- KMO(mat)
overall_msa <- round(kmo_res$MSA, 3)
msa_items <- data.frame(variable = names(kmo_res$MSAi),
                        MSA = round(kmo_res$MSAi, 3)) %>%
       arrange(MSA)
cat("Overall MSA =", overall_msa, "\n\n")
## Overall MSA = 0.669
print(msa_items)
##          variable   MSA
## X7_resid X7_resid 0.104
## X10           X10 0.306
## X8             X8 0.401
## X9             X9 0.457
## X6             X6 0.514
## X1             X1 0.547
## X2             X2 0.728
## X4             X4 0.814
## X5             X5 0.859
## X3             X3 0.870
print(kmo_res)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = mat)
## Overall MSA =  0.67
## MSA for each item = 
##       X1       X2       X3       X4       X5       X6 X7_resid       X8 
##     0.55     0.73     0.87     0.81     0.86     0.51     0.10     0.40 
##       X9      X10 
##     0.46     0.31

Hasil uji KMO menunjukkan nilai Overall MSA sebesar 0,669 (≈0,67) yang tergolong cukup (mediocre) sehingga data masih layak untuk analisis faktor, namun tidak optimal. Secara per item, variabel X3 (0,870), X5 (0,859), X4 (0,814), dan X2 (0,728) memiliki MSA baik sehingga sangat mendukung pembentukan faktor. Sebaliknya, X7_resid (0,104), X10 (0,306), X8 (0,401), dan X9 (0,457) memiliki MSA rendah (<0,5) yang menunjukkan kontribusinya terhadap struktur faktor kurang memadai, sementara X6 (0,514) dan X1 (0,547) berada pada batas minimal kelayakan. Dengan demikian, walaupun model secara keseluruhan masih dapat dianalisis, namun masih terdapat beberapa variabel dengan MSA rendah yang sebaiknya dipertimbangkan untuk dieliminasi agar kualitas faktor meningkat.

Uji Bartlett

bartlett_res <- cortest.bartlett(mat, n = nrow(df_log))
print(bartlett_res) # chi-square dan p-value
## $chisq
## [1] 129842.7
## 
## $p.value
## [1] 0
## 
## $df
## [1] 45

Hasil uji Bartlett ini menunjukkan nilai chi-square sebesar 129842,7 dengan derajat kebebasan 45 dan p-value = 0 (<0,05), sehingga H0 yang menyatakan matriks korelasi adalah matriks identitas ditolak. Artinya, terdapat korelasi yang signifikan antar variabel dan data memenuhi syarat untuk dilakukan analisis faktor atau PCA karena struktur korelasinya yang tidak bersifat acak.

PCA (Principal Component Analysis)

Pada studi kasus ini, kami menerapkan kombinasi dua pendekatan PCA, yaitu manual (prcomp + eigen) dan fungsi principal() untuk rotasi dan interpretasi faktor.

Menjalankan PCA & Eigenvalue

vars_pca <- c("X1","X2","X3","X4","X5","X6","X7_resid","X8","X9","X10")
X <- df_log[vars_pca]
scale_data <- scale(X)
pca <- prcomp(scale_data, center = FALSE, scale. = FALSE)
explained <- (pca$sdev^2)/sum(pca$sdev^2)
eig_table <- data.frame(PC = paste0("PC",1:length(explained)),
                               eigenvalue = round(pca$sdev^2,3),
                               propvar = round(explained,3),
                               cumvar = round(cumsum(explained),3))
print(eig_table)
##      PC eigenvalue propvar cumvar
## 1   PC1      3.850   0.385  0.385
## 2   PC2      1.163   0.116  0.501
## 3   PC3      1.013   0.101  0.603
## 4   PC4      1.001   0.100  0.703
## 5   PC5      0.987   0.099  0.801
## 6   PC6      0.861   0.086  0.887
## 7   PC7      0.400   0.040  0.927
## 8   PC8      0.365   0.037  0.964
## 9   PC9      0.275   0.028  0.991
## 10 PC10      0.085   0.009  1.000

Pada eigenvalue di atas, PC1 memiliki eigenvalue 3,850 dan menjelaskan 38,5% variasi data. PC2 dan PC3 masing-masing menjelaskan 11,6% dan 10,1%. Berdasarkan kriteria Kaiser (eigenvalue > 1), terdapat 4 komponen yang layak dipertahankan (PC1–PC4). Namun, 3 komponen pertama sudah menjelaskan 60,3% variasi kumulatif sehingga secara praktis 3 komponen cukup merepresentasikan struktur utama data.

Scree Plot

fviz_eig(pca, addlabels = TRUE)
## Warning in geom_bar(stat = "identity", fill = barfill, color = barcolor, :
## Ignoring empty aesthetic: `width`.

Berdasarkan scree plot, terlihat bahwa penurunan yang sangat tajam dari PC1 (38,5%) ke PC2 (11,6%), kemudian penurunan mulai melandai pada PC3 (10,1%) dan seterusnya relatif stabil di sekitar 10% hingga di bawahnya. Pola ini menunjukkan adanya titik “elbow” sekitar komponen ke-3, yang menandakan bahwa tambahan komponen setelah PC3 hanya memberikan peningkatan variasi yang relatif kecil. Dengan demikian, secara visual scree plot mendukung pemilihan 3 komponen utama sebagai representasi yang cukup untuk menjelaskan struktur variansi data.

Menyimpan Skor PCA

scores_pca <- as.data.frame(pca$x)
df_log <- cbind(df_log, scores_pca)

Setiap baris data memiliki nilai PC1–PC10 yang merupakan representasi baru dalam ruang komponen utama, sehingga dapat digunakan untuk analisis lanjutan.

Rotasi Varimax

k <- 3
Lk <- pca$rotation[, 1:k]
vm <- varimax(Lk)
rot_loadings_varimax <- as.matrix(vm$loadings)
scores_varimax <- as.matrix(scale_data) %*% rot_loadings_varimax
colnames(scores_varimax) <- paste0("VM_PC", 1:k)
df_log <- cbind(df_log, as.data.frame(scores_varimax))
L_vm <- as.matrix(rot_loadings_varimax)
cat("\nVarimax loadings (full precision):\n")
## 
## Varimax loadings (full precision):
print(formatC(L_vm, digits = 4, format = "f"), quote = FALSE)
##          PC1     PC2     PC3    
## X1       0.0000  -0.0045 0.7221 
## X2       -0.4184 0.1022  -0.0186
## X3       -0.4382 -0.0054 0.0015 
## X4       -0.4561 0.0124  0.0018 
## X5       -0.4278 0.0218  0.0084 
## X6       -0.4326 -0.3377 0.0456 
## X7_resid -0.0693 -0.6963 -0.0438
## X8       -0.2234 0.6232  -0.0523
## X9       0.0048  0.0378  0.6629 
## X10      -0.0051 -0.0148 0.1789

Rotasi Varimax menunjukkan bahwa PC1 didominasi oleh X2–X6 (loading sekitar −0,41 sampai −0,46), PC2 oleh X7_resid (−0,6963) dan X8 (0,6232), serta PC3 oleh X1 (0,7221) dan X9 (0,6629). Tidak terdapat cross-loading signifikan (≥0,4 pada lebih dari satu komponen), sehingga struktur faktor tergolong bersih. Variabel X10 memiliki loading rendah pada semua komponen, sehingga kontribusinya dalam model relatif lemah.Secara umum, tidak terlihat cross-loading yang berarti (≥0,4 pada lebih dari satu komponen), sehingga struktur hasil rotasi Varimax dapat dikatakan cukup bersih dan interpretatif.

Rotasi Oblimin

pc_oblimin <- principal(scale_data, nfactors = k, rotate = "oblimin", scores = TRUE)
## Loading required namespace: GPArotation
rot_loadings_oblimin <- as.matrix(pc_oblimin$loadings)
scores_oblimin <- as.data.frame(pc_oblimin$scores)
colnames(scores_oblimin) <- paste0("OB_PC", 1:k)
df_log <- cbind(df_log, scores_oblimin)
L_ob <- as.matrix(rot_loadings_oblimin)
cat("\nOblimin loadings (full precision):\n")
## 
## Oblimin loadings (full precision):
print(formatC(L_ob, digits = 4, format = "f"), quote = FALSE)
##          TC1    TC2     TC3    
## X1       0.0241 -0.0003 0.7273 
## X2       0.7970 0.1653  -0.0220
## X3       0.8443 0.0486  -0.0040
## X4       0.8773 0.0706  -0.0036
## X5       0.8223 0.0774  0.0036 
## X6       0.8608 -0.3206 0.0344 
## X7_resid 0.1865 -0.7643 -0.0577
## X8       0.3797 0.7189  -0.0441
## X9       0.0095 0.0456  0.6685 
## X10      0.0169 -0.0147 0.1799

Hasil rotasi Oblimin di atas menunjukkan struktur faktor yang lebih kuat dan terkelompok jelas. TC1 didominasi oleh X2–X6 dengan loading sangat tinggi (0,79–0,88), sehingga merepresentasikan satu dimensi utama yang kuat. TC2 terutama dibentuk oleh X7_resid (−0,7643) dan X8 (0,7189), sementara TC3 didominasi oleh X1 (0,7273) dan X9 (0,6685). Tidak terdapat cross-loading yang signifikan (≥0,4 pada lebih dari satu faktor), sehingga struktur tetap bersih meskipun rotasi bersifat oblique. Variabel X10 kembali menunjukkan loading rendah pada semua faktor, sehingga kontribusinya relatif lemah dalam model.

Evaluasi Struktur Loading & Communalities

Communalities (k = 3 Komponen)

L_vm <- as.matrix(rot_loadings_varimax)
L_ob <- as.matrix(rot_loadings_oblimin)
comm_varimax_full <- rowSums(L_vm^2)
comm_oblimin_full <- rowSums(L_ob^2)
comm_table_full <- data.frame(variable = rownames(L_vm),
                                 comm_varimax = round(comm_varimax_full, 3),
                                 comm_oblimin = round(comm_oblimin_full, 3))
cat("\nCommunalities (k components):\n")
## 
## Communalities (k components):
print(comm_table_full, row.names = FALSE)
##  variable comm_varimax comm_oblimin
##        X1        0.521        0.530
##        X2        0.186        0.663
##        X3        0.192        0.715
##        X4        0.208        0.775
##        X5        0.184        0.682
##        X6        0.303        0.845
##  X7_resid        0.492        0.622
##        X8        0.441        0.663
##        X9        0.441        0.449
##       X10        0.032        0.033

Untuk k = 3 komponen, rotasi oblimin menghasilkan nilai komunalitas yang lebih tinggi dibandingkan varimax, sehingga lebih mampu menjelaskan varians sebagian besar variabel (terutama X2–X6). Pada varimax, banyak variabel memiliki komunalitas rendah, sedangkan X10 sangat rendah pada kedua metode, menandakan kurang terwakili oleh model. Secara keseluruhan, oblimin lebih optimal dalam merepresentasikan struktur data.

Maximum Absolute Loading per Variable

max_abs_vm <- apply(abs(L_vm), 1, max)
max_abs_ob <- apply(abs(L_ob), 1, max)
max_table <- data.frame(variable = rownames(L_vm),
                           max_abs_varimax = round(max_abs_vm, 3),
                           max_abs_oblimin = round(max_abs_ob, 3))
cat("\nMax absolute loading per variable:\n")
## 
## Max absolute loading per variable:
print(max_table, row.names = FALSE)
##  variable max_abs_varimax max_abs_oblimin
##        X1           0.722           0.727
##        X2           0.418           0.797
##        X3           0.438           0.844
##        X4           0.456           0.877
##        X5           0.428           0.822
##        X6           0.433           0.861
##  X7_resid           0.696           0.764
##        X8           0.623           0.719
##        X9           0.663           0.668
##       X10           0.179           0.180

Tabel maximum absolute loading menunjukkan kekuatan hubungan terbesar setiap variabel dengan salah satu dari tiga komponen. Terlihat pada rotasi oblimin, hampir semua variabel (terutama X2–X6) memiliki loading maksimum yang jauh lebih tinggi (≥ 0,79 hingga 0,88) dibandingkan varimax yang sebagian besar berada di kisaran 0,41–0,46. Hal ini menandakan bahwa pada oblimin, struktur faktor terlihat lebih tegas dan variabel lebih kuat terasosiasi dengan komponennya masing-masing. Sementara itu, X10 memiliki loading sangat rendah pada kedua metode (~0,18), sehingga kontribusinya terhadap komponen yang terbentuk relatif lemah.

Matriks Loading Numerik

L_vm_num <- as.matrix(rot_loadings_varimax)
L_ob_num <- as.matrix(rot_loadings_oblimin)
print(as.data.frame(formatC(L_vm_num, digits = 3, format = "f")), quote = FALSE)
##             PC1    PC2    PC3
## X1        0.000 -0.004  0.722
## X2       -0.418  0.102 -0.019
## X3       -0.438 -0.005  0.001
## X4       -0.456  0.012  0.002
## X5       -0.428  0.022  0.008
## X6       -0.433 -0.338  0.046
## X7_resid -0.069 -0.696 -0.044
## X8       -0.223  0.623 -0.052
## X9        0.005  0.038  0.663
## X10      -0.005 -0.015  0.179

Matriks loadings menunjukkan bahwa setiap variabel cenderung memiliki satu loading dominan pada salah satu dari tiga komponen. PC1 terutama memuat X2–X6 dengan loading negatif moderat (sekitar −0,42 hingga −0,46), PC2 didominasi oleh X7_resid (−0,696) dan X8 (0,623), sedangkan PC3 kuat pada X1 (0,722) dan X9 (0,663). Sebagian besar variabel memiliki loading kecil pada komponen lainnya, menandakan struktur yang cukup bersih (minim cross-loading). Variabel X10 memiliki loading rendah di semua komponen (maksimum 0,179), sehingga kontribusinya terhadap struktur tiga komponen relatif lemah.

Annotasi Loading Signifikan (|loading| ≥ 0.4)

threshold <- 0.4
L_vm_annot <- matrix(sprintf("%.6f", L_vm_num), nrow = nrow(L_vm_num), dimnames = dimnames(L_vm_num))
L_ob_annot <- matrix(sprintf("%.6f", L_ob_num), nrow = nrow(L_ob_num), dimnames = dimnames(L_ob_num))
L_vm_annot <- ifelse(abs(L_vm_num) >= threshold, paste0(L_vm_annot, " *"), L_vm_annot)
L_ob_annot <- ifelse(abs(L_ob_num) >= threshold, paste0(L_ob_annot, " *"), L_ob_annot)
cat("\nVarimax loadings annotated (\"*\" = |loading| >= 0.4):\n")
## 
## Varimax loadings annotated ("*" = |loading| >= 0.4):
print(as.data.frame(L_vm_annot), quote = FALSE)
##                  PC1         PC2        PC3
## X1          0.000000   -0.004452 0.722111 *
## X2       -0.418370 *    0.102175  -0.018563
## X3       -0.438225 *   -0.005372   0.001461
## X4       -0.456110 *    0.012440   0.001804
## X5       -0.427798 *    0.021755   0.008447
## X6       -0.432559 *   -0.337674   0.045559
## X7_resid   -0.069340 -0.696303 *  -0.043837
## X8         -0.223413  0.623196 *  -0.052341
## X9          0.004827    0.037787 0.662863 *
## X10        -0.005114   -0.014844   0.178914
cat("\nOblimin loadings annotated (\"*\" = |loading| >= 0.4):\n")
## 
## Oblimin loadings annotated ("*" = |loading| >= 0.4):
print(as.data.frame(L_ob_annot), quote = FALSE)
##                 TC1         TC2        TC3
## X1         0.024063   -0.000269 0.727339 *
## X2       0.796975 *    0.165297  -0.022002
## X3       0.844277 *    0.048563  -0.004033
## X4       0.877333 *    0.070554  -0.003584
## X5       0.822307 *    0.077411   0.003626
## X6       0.860814 *   -0.320577   0.034412
## X7_resid   0.186549 -0.764273 *  -0.057687
## X8         0.379712  0.718950 *  -0.044140
## X9         0.009519    0.045615 0.668483 *
## X10        0.016884   -0.014677   0.179896

Berdasarkan loading yang diberi tanda (*) untuk |loading| ≥ 0,4, terlihat bahwa pada kedua rotasi setiap variabel (kecuali X10) memiliki satu loading dominan sehingga tidak terjadi cross-loading. Pada varimax, PC1 memuat kuat X2–X6, PC2 memuat X7_resid dan X8, serta PC3 memuat X1 dan X9. Pada oblimin, pola yang terbentuk serupa tetapi dengan nilai loading yang lebih tinggi dan lebih tegas, terutama pada X2–X6 di TC1 serta X7_resid dan X8 di TC2. X10 tidak mencapai ambang 0,4 pada semua komponen, sehingga kontribusinya terhadap struktur faktor relatif lemah.

PCA Manual (Eigen Decomposition)

Matriks Kovarian Data Terskalakan

R <- cov(scale_data)

Eigen Decomposition Manual

pc_manual <- eigen(R)
evalues <- pc_manual$values
evectors <- pc_manual$vectors
sumvar <- sum(evalues)
propvar_manual <- evalues / sumvar
cumvar_manual <- cumsum(propvar_manual)
cumvar_df <- data.frame(eigen_value = round(evalues, 6),
                           propvar = round(propvar_manual, 6),
                           cumvar = round(cumvar_manual, 6))
row.names(cumvar_df) <- paste0("PC", 1:ncol(scale_data))
print(cumvar_df)
##      eigen_value  propvar   cumvar
## PC1     3.849792 0.384979 0.384979
## PC2     1.163277 0.116328 0.501307
## PC3     1.013008 0.101301 0.602608
## PC4     1.001113 0.100111 0.702719
## PC5     0.986530 0.098653 0.801372
## PC6     0.861255 0.086125 0.887498
## PC7     0.399667 0.039967 0.927464
## PC8     0.365040 0.036504 0.963968
## PC9     0.275262 0.027526 0.991494
## PC10    0.085056 0.008506 1.000000

Hasil eigen decomposition menunjukkan bahwa PC1 memiliki eigenvalue terbesar (3,85) dan menjelaskan sekitar 38,5% variasi data, diikuti PC2 (11,6%) dan PC3 (10,1%), sehingga tiga komponen pertama secara kumulatif menjelaskan sekitar 60,3% total varians. Hingga PC4, varians kumulatif mencapai 70,3%. Berdasarkan kriteria Kaiser (eigenvalue > 1), terdapat empat komponen yang layak dipertahankan (PC1–PC4), meskipun PC3 dan PC4 nilainya sangat dekat dengan 1. Komponen setelah PC4 memiliki kontribusi varians relatif kecil, sehingga perannya dalam menjelaskan struktur data semakin terbatas.

Validasi PCA Manual vs prcomp()

Bandingkan Eigenvector

rot_prcomp <- pca$rotation
cor_mat <- abs(cor(evectors, rot_prcomp))
print(round(cor_mat, 3))
##         PC1   PC2   PC3   PC4   PC5   PC6   PC7 PC8   PC9  PC10
##  [1,] 1.000 0.202 0.649 0.249 0.115 0.370 0.007   0 0.004 0.138
##  [2,] 0.202 1.000 0.096 0.037 0.017 0.055 0.001   0 0.001 0.020
##  [3,] 0.649 0.096 1.000 0.118 0.055 0.176 0.003   0 0.002 0.066
##  [4,] 0.249 0.037 0.118 1.000 0.021 0.068 0.001   0 0.001 0.025
##  [5,] 0.115 0.017 0.055 0.021 1.000 0.031 0.001   0 0.000 0.012
##  [6,] 0.370 0.055 0.176 0.068 0.031 1.000 0.002   0 0.001 0.037
##  [7,] 0.007 0.001 0.003 0.001 0.001 0.002 1.000   0 0.000 0.001
##  [8,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000   1 0.000 0.000
##  [9,] 0.004 0.001 0.002 0.001 0.000 0.001 0.000   0 1.000 0.000
## [10,] 0.138 0.020 0.066 0.025 0.012 0.037 0.001   0 0.000 1.000

Matriks korelasi absolut antara eigenvector manual dan hasil prcomp menunjukkan nilai 1,000 pada diagonal utama untuk setiap pasangan komponen yang bersesuaian (misalnya PC1 dengan PC1, PC2 dengan PC2, dan seterusnya), yang berarti arah vektor identik (atau sama hingga perbedaan tanda). Nilai di luar diagonal relatif kecil, menandakan tidak ada pertukaran atau pencampuran komponen. Dengan demikian, hasil perhitungan eigen secara manual konsisten dan setara dengan hasil PCA menggunakan prcomp.

Hitung Skor PCA Manual

scores_manual <- as.matrix(scale_data) %*% evectors
head(round(scores_manual[,1:3], 4))
##         [,1]    [,2]    [,3]
## [1,]  0.1461 -0.4258 -0.5120
## [2,] -1.8925  1.6221  1.9102
## [3,]  1.8864 -0.0833  0.2806
## [4,]  2.9448  1.5138  0.0811
## [5,]  0.4821 -1.6749  0.2658
## [6,]  0.8949  0.8818 -0.7168
head(round(pca$x[,1:3], 4))
##          PC1     PC2     PC3
## [1,]  0.1461 -0.4258 -0.5120
## [2,] -1.8925  1.6221  1.9102
## [3,]  1.8864 -0.0833  0.2806
## [4,]  2.9448  1.5138  0.0811
## [5,]  0.4821 -1.6749  0.2658
## [6,]  0.8949  0.8818 -0.7168

Cek Korelasi Skor

print(round(cor(scores_manual[,1:3], pca$x[,1:3]), 3))
##      PC1 PC2 PC3
## [1,]   1   0   0
## [2,]   0   1   0
## [3,]   0   0   1

Evaluasi Ulang Communalities

print(comm_table_full)
##          variable comm_varimax comm_oblimin
## X1             X1        0.521        0.530
## X2             X2        0.186        0.663
## X3             X3        0.192        0.715
## X4             X4        0.208        0.775
## X5             X5        0.184        0.682
## X6             X6        0.303        0.845
## X7_resid X7_resid        0.492        0.622
## X8             X8        0.441        0.663
## X9             X9        0.441        0.449
## X10           X10        0.032        0.033

Cross-Loadings Check

threshold <- 0.4
cross_vm <- apply(abs(L_vm_num) >= threshold, 1, sum)
cross_ob <- apply(abs(L_ob_num) >= threshold, 1, sum)
data.frame(variable = rownames(L_vm_num), cross_varimax = cross_vm, cross_oblimin = cross_ob)
##          variable cross_varimax cross_oblimin
## X1             X1             1             1
## X2             X2             1             1
## X3             X3             1             1
## X4             X4             1             1
## X5             X5             1             1
## X6             X6             1             1
## X7_resid X7_resid             1             1
## X8             X8             1             1
## X9             X9             1             1
## X10           X10             0             0

Hasil cross-loading menunjukkan bahwa hampir semua variabel (X1–X9) memiliki tepat satu loading yang ≥ 0,4 pada masing-masing rotasi, baik varimax maupun oblimin, sehingga tidak terjadi cross-loading (tidak ada variabel yang kuat pada lebih dari satu komponen). Ini menandakan struktur faktor yang bersih dan mudah diinterpretasikan. Sementara itu, X10 memiliki nilai 0 pada kedua metode, artinya tidak ada loading yang mencapai ambang 0,4 sehingga variabel tersebut tidak berkontribusi kuat pada komponen mana pun.

Korelasi Antar Faktor (Oblimin)

round(cor(scores_oblimin, use = "pairwise.complete.obs"), 3)
##        OB_PC1 OB_PC2 OB_PC3
## OB_PC1  1.000  0.168 -0.036
## OB_PC2  0.168  1.000 -0.027
## OB_PC3 -0.036 -0.027  1.000

Hasil korelasi antar faktor dengan rotasi oblimin menunjukkan bahwa hubungan antar faktor relatif rendah. Korelasi antara Faktor 1 dan Faktor 2 sebesar 0,168 mengindikasikan hubungan yang lemah, sedangkan korelasi antara Faktor 1 dan Faktor 3 (-0,036) serta Faktor 2 dan Faktor 3 (-0,027) mendekati nol. Hal ini menunjukkan bahwa ketiga faktor yang terbentuk cenderung saling independen meskipun menggunakan rotasi oblique (oblimin), sehingga struktur faktor dapat dikatakan stabil dan tidak terjadi tumpang tindih yang kuat antar dimensi.

Visualisasi

Biplot

p <- fviz_pca_biplot(
         pca,
         geom.ind = "point",    # gambar titik individu
         col.ind = "orange",    # warna titik
         pointsize = 2.0,       # ukuran titik
         alpha.ind = 0.35,      # transparansi titik (0..1)
         repel = TRUE,          # hindari label saling tumpuk
         label = "var"          # tampilkan panah/label variabel
     )

p

Visualisasi PCA biplot menunjukkan proyeksi seluruh observasi sebagai titik oranye dan kontribusi masing-masing variabel sebagai panah pada dua komponen utama, yaitu Dim1 (38,5%) dan Dim2 (11,6%). Sebaran titik yang lebih melebar secara horizontal menandakan bahwa variasi data lebih dominan dijelaskan oleh Dim1. Panah-panah variabel yang mengarah ke kanan dan relatif sejajar menunjukkan adanya korelasi positif antar variabel tersebut serta kontribusi kuat terhadap Komponen Utama 1, sementara variabel dengan arah berbeda (ke atas atau ke bawah) lebih berkontribusi pada Dim2 atau memiliki pola hubungan yang berbeda. Panjang panah mencerminkan besar kontribusi variabel dalam membentuk struktur data, sehingga variabel dengan panah lebih panjang memiliki pengaruh yang lebih dominan dalam menjelaskan dinamika performa akun Instagram berdasarkan metrik engagement dan exposure.

Plot Circle

contrib_circle <- fviz_pca_var(pca, col.var = "contrib",
                                                                   gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                                                                   repel = TRUE) + ggtitle("Kontribusi Variabel")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(contrib_circle)

Berdasarkan grafik correlation circle hasil fviz_pca_var(), terlihat bahwa dua komponen utama pertama mampu menjelaskan sebesar 50,1% variasi data (Dim1 = 38,5% dan Dim2 = 11,6%). Variabel dengan kontribusi paling besar terhadap pembentukan komponen utama adalah X8 dan X7_resid, diikuti oleh X6 dan X3, yang ditunjukkan oleh panah lebih panjang dan warna lebih pekat. Variabel X2, X3, X4, dan X5 memiliki arah yang relatif searah sehingga mengindikasikan adanya korelasi positif antarvariabel tersebut, sedangkan X7_resid memiliki arah yang berbeda (negatif pada Dim2) sehingga menunjukkan pola hubungan yang berlawanan dengan beberapa variabel lain. Secara umum, variabel yang berada dekat dengan lingkaran menunjukkan representasi yang baik pada dua komponen utama pertama.

Grafik Kontribusi

contrib_v_PC1 <- fviz_contrib(pca, choice = "var", axes = 1, top = 5) + ggtitle("PC1")
contrib_v_PC2 <- fviz_contrib(pca, choice = "var", axes = 2, top = 5) + ggtitle("PC2")
contrib_v_PC3 <- fviz_contrib(pca, choice = "var", axes = 3, top = 5) + ggtitle("PC3")
plot(contrib_v_PC1)

plot(contrib_v_PC2)

plot(contrib_v_PC3)

Berdasarkan grafik kontribusi variabel terhadap tiga komponen utama, terlihat bahwa PC1 terutama dibentuk oleh variabel X4, X3, X2, X5, dan X6 dengan kontribusi relatif merata (sekitar 15–21%) dan seluruhnya berada di atas garis rata-rata (10%), sehingga PC1 merepresentasikan kombinasi kuat dari kelima variabel tersebut. Pada PC2, kontribusi paling dominan berasal dari X7_resid (≈49%) dan X8 (≈35%), sedangkan X6 berkontribusi sedang dan variabel lain sangat kecil, sehingga PC2 terutama mencerminkan variasi yang dipengaruhi kuat oleh X7_resid. Sementara itu, PC3 sangat didominasi oleh X1 (≈52%) dan X9 (≈44%), dengan kontribusi variabel lain sangat minimal, yang menunjukkan bahwa komponen ketiga terutama menangkap pola variasi spesifik dari dua variabel tersebut. Secara keseluruhan, setiap komponen utama memiliki kelompok variabel dominan yang berbeda, sehingga masing-masing merepresentasikan dimensi variasi yang berbeda dalam data.

FA (Factor Analysis)

Uji Kelayakan Awal (KMO & Bartlett)

kmo_res <- KMO(scale_data)
print(kmo_res)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = scale_data)
## Overall MSA =  0.67
## MSA for each item = 
##       X1       X2       X3       X4       X5       X6 X7_resid       X8 
##     0.55     0.73     0.87     0.81     0.86     0.51     0.10     0.40 
##       X9      X10 
##     0.46     0.31
bartlett_res <- cortest.bartlett(cor(scale_data), n = nrow(scale_data))
print(bartlett_res)
## $chisq
## [1] 129842.7
## 
## $p.value
## [1] 0
## 
## $df
## [1] 45

Membuat Matriks Korelasi & Eigen Decomposition (10 Variabel Awal)

head(scale_data, 6)
##               X1          X2          X3         X4         X5         X6
## [1,] -1.46633158  0.07206262 -0.06765457 -0.2779544  0.2353923 -0.1800679
## [2,]  2.14893794  0.79593216  0.59063940  0.8142224  0.8611206  0.5959253
## [3,] -0.07873629 -0.38504092 -0.82044727 -1.7746696 -0.1516591 -1.5659035
## [4,]  0.09121191 -0.57819752 -2.01359546 -0.2779544 -3.0421829 -0.7627444
## [5,]  1.04030694 -0.12666066  0.37270093 -0.5885508 -0.1926380  0.1229421
## [6,] -1.46633158 -0.15504405 -0.50801007 -0.7853946 -0.1124246 -0.8704109
##        X7_resid         X8          X9        X10
## [1,]  0.6521639 -0.1414844  0.90433325 -0.1835984
## [2,] -1.7607805  1.0238640  0.79138293 -1.0277124
## [3,]  1.6104745  0.4743697  0.72986344  0.1619995
## [4,] -1.8173304 -0.4655013  0.08165148 -0.1835984
## [5,]  1.5652196 -0.8436764 -0.89013681  0.4711474
## [6,] -0.2804005  0.4114918  0.91000599 -1.0277124
R <- cor(scale_data)
eig <- eigen(R)
cat("Eigen values:\n")
## Eigen values:
print(round(eig$values, 4))
##  [1] 3.8498 1.1633 1.0130 1.0011 0.9865 0.8613 0.3997 0.3650 0.2753 0.0851
cat("\nEigen vectors:\n")
## 
## Eigen vectors:
print(round(eig$vectors, 4))
##          [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
##  [1,]  0.0052  0.0137  0.7220  0.0897  0.6828 -0.0656 -0.0004 -0.0009 -0.0050
##  [2,] -0.4281  0.0471 -0.0170  0.0004  0.0104 -0.1364  0.7539  0.2587  0.2304
##  [3,] -0.4339 -0.0615  0.0058 -0.0015  0.0034  0.0221 -0.1947 -0.7356  0.4488
##  [4,] -0.4539 -0.0462  0.0059 -0.0076 -0.0060  0.0028 -0.0804 -0.1439 -0.8237
##  [5,] -0.4270 -0.0331  0.0120  0.0060 -0.0073 -0.0045 -0.5991  0.6044  0.2510
##  [6,] -0.3853 -0.3891  0.0578 -0.0197 -0.0033  0.4809  0.1640  0.0683 -0.0259
##  [7,]  0.0205 -0.7003 -0.0264 -0.0124 -0.0245 -0.7025 -0.0152 -0.0053 -0.0135
##  [8,] -0.3020  0.5879 -0.0652  0.0192  0.0070 -0.4992 -0.0314 -0.0358 -0.0560
##  [9,]  0.0042  0.0547  0.6617 -0.3475 -0.6600 -0.0509  0.0121 -0.0038  0.0075
## [10,] -0.0020 -0.0109  0.1793  0.9329 -0.3120 -0.0014  0.0110 -0.0056 -0.0040
##         [,10]
##  [1,] -0.0018
##  [2,]  0.3273
##  [3,]  0.1645
##  [4,]  0.2934
##  [5,]  0.1708
##  [6,] -0.6579
##  [7,] -0.1169
##  [8,] -0.5514
##  [9,]  0.0022
## [10,] -0.0017

Hasil eigen decomposition menunjukkan terdapat empat eigenvalue > 1 (3,8498; 1,1633; 1,0130; 1,0011), sehingga berdasarkan kriteria Kaiser dapat dipertahankan empat komponen utama. PC1 memiliki eigenvalue terbesar (3,8498), artinya menjelaskan variasi paling dominan dalam data, sedangkan komponen berikutnya menjelaskan variasi tambahan yang lebih kecil. Nilai pada eigenvector menunjukkan bobot masing-masing variabel dalam membentuk setiap komponen, di mana semakin besar nilai absolutnya, semakin besar kontribusinya terhadap komponen tersebut.

Refinement Tahap 1 (Mengurangi Variabel)

vars_refined <- c("X2","X3","X4","X5","X6","X7_resid","X8")
data_refined <- scale(df_log[vars_refined])
KMO(data_refined)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data_refined)
## Overall MSA =  0.67
## MSA for each item = 
##       X2       X3       X4       X5       X6 X7_resid       X8 
##     0.73     0.87     0.81     0.86     0.51     0.10     0.40
cortest.bartlett(cor(data_refined), n = nrow(data_refined))
## $chisq
## [1] 129814.1
## 
## $p.value
## [1] 0
## 
## $df
## [1] 21

Pada tahap refinement pertama, nilai KMO keseluruhan sebesar 0,67 menunjukkan bahwa data sudah cukup layak untuk analisis faktor (kategori sedang), namun nilai MSA per variabel menunjukkan bahwa X7_resid (0,10) dan X8 (0.40) memiliki kecukupan sampel yang rendah sehingga kurang mendukung pembentukan faktor. Sementara itu, variabel lain seperti X2, X3, X4, dan X5 memiliki MSA yang baik (>0,7). Hasil uji Bartlett menghasilkan p-value = 0 (signifikan), yang berarti matriks korelasi tidak berbentuk identity matrix, sehingga antar variabel terdapat korelasi yang cukup untuk dilakukan analisis faktor.

Refinement Tahap 2 (Variabel Final)

vars_final <- c("X2","X3","X4","X5","X6")
data_final <- scale(df_log[vars_final])
KMO(data_final)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data_final)
## Overall MSA =  0.89
## MSA for each item = 
##   X2   X3   X4   X5   X6 
## 0.91 0.89 0.86 0.90 0.91
cortest.bartlett(cor(data_final), n = nrow(data_final))
## $chisq
## [1] 86729.72
## 
## $p.value
## [1] 0
## 
## $df
## [1] 10

Pada tahap refinement kedua, nilai KMO keseluruhan meningkat menjadi 0,89, yang menunjukkan tingkat kelayakan yang sangat baik untuk analisis faktor. Seluruh variabel (X2–X6) juga memiliki nilai MSA tinggi (0,86–0,91), menandakan bahwa masing-masing variabel sudah memadai dalam membentuk struktur faktor. Hasil uji Bartlett tetap signifikan (p-value = 0), sehingga terdapat korelasi yang cukup antar variabel dan analisis faktor sangat layak untuk dilanjutkan pada model final ini.

Menentukan Jumlah Faktor (Parallel Analysis)

R <- cor(data_final)
eig <- eigen(R)
print(round(eig$values,4))
## [1] 3.5701 0.4131 0.3879 0.3577 0.2713
fa.parallel(data_final, fa="fa", fm="ml")

## Parallel analysis suggests that the number of factors =  1  and the number of components =  NA

Estimasi Factor Analysis

2 Faktor Tanpa Rotasi

fa_final <- fa(
       r = data_final,
       nfactors = 2,
       rotate = "none",
       fm = "ml"
   )
print(fa_final)
## Factor Analysis using method =  ml
## Call: fa(r = data_final, nfactors = 2, rotate = "none", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     ML1   ML2   h2    u2 com
## X2 0.67  0.36 0.58 0.423 1.5
## X3 0.72  0.38 0.66 0.343 1.5
## X4 1.00 -0.01 1.00 0.005 1.0
## X5 0.69  0.37 0.61 0.386 1.5
## X6 0.68  0.38 0.61 0.388 1.6
## 
##                        ML1  ML2
## SS loadings           2.90 0.56
## Proportion Var        0.58 0.11
## Cumulative Var        0.58 0.69
## Proportion Explained  0.84 0.16
## Cumulative Proportion 0.84 1.00
## 
## Mean item complexity =  1.4
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  10  with the objective function =  2.89 with Chi Square =  86729.72
## df of  the model are 1  and the objective function was  0 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0.01 
## 
## The harmonic n.obs is  29999 with the empirical chi square  1.28  with prob <  0.26 
## The total n.obs was  29999  with Likelihood Chi Square =  8.43  with prob <  0.0037 
## 
## Tucker Lewis Index of factoring reliability =  0.999
## RMSEA index =  0.016  and the 90 % confidence intervals are  0.007 0.026
## BIC =  -1.88
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1  ML2
## Correlation of (regression) scores with factors   1.00 0.77
## Multiple R square of scores with factors          1.00 0.60
## Minimum correlation of possible factor scores     0.99 0.20

Estimasi Ulang

fa_2_none <- fa(
       r = data_final,
       nfactors = 2,
       rotate = "none",
       fm = "ml"
   )
print(fa_2_none)
## Factor Analysis using method =  ml
## Call: fa(r = data_final, nfactors = 2, rotate = "none", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     ML1   ML2   h2    u2 com
## X2 0.67  0.36 0.58 0.423 1.5
## X3 0.72  0.38 0.66 0.343 1.5
## X4 1.00 -0.01 1.00 0.005 1.0
## X5 0.69  0.37 0.61 0.386 1.5
## X6 0.68  0.38 0.61 0.388 1.6
## 
##                        ML1  ML2
## SS loadings           2.90 0.56
## Proportion Var        0.58 0.11
## Cumulative Var        0.58 0.69
## Proportion Explained  0.84 0.16
## Cumulative Proportion 0.84 1.00
## 
## Mean item complexity =  1.4
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  10  with the objective function =  2.89 with Chi Square =  86729.72
## df of  the model are 1  and the objective function was  0 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0.01 
## 
## The harmonic n.obs is  29999 with the empirical chi square  1.28  with prob <  0.26 
## The total n.obs was  29999  with Likelihood Chi Square =  8.43  with prob <  0.0037 
## 
## Tucker Lewis Index of factoring reliability =  0.999
## RMSEA index =  0.016  and the 90 % confidence intervals are  0.007 0.026
## BIC =  -1.88
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1  ML2
## Correlation of (regression) scores with factors   1.00 0.77
## Multiple R square of scores with factors          1.00 0.60
## Minimum correlation of possible factor scores     0.99 0.20

Rotasi Varimax

fa_2_varimax <- fa(
       r = data_final,
       nfactors = 2,
       rotate = "varimax",
       fm = "ml"
   )
print(fa_2_varimax)
## Factor Analysis using method =  ml
## Call: fa(r = data_final, nfactors = 2, rotate = "varimax", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     ML1  ML2   h2    u2 com
## X2 0.53 0.54 0.58 0.423 2.0
## X3 0.57 0.58 0.66 0.343 2.0
## X4 0.95 0.29 1.00 0.005 1.2
## X5 0.54 0.56 0.61 0.386 2.0
## X6 0.54 0.57 0.61 0.388 2.0
## 
##                        ML1  ML2
## SS loadings           2.10 1.36
## Proportion Var        0.42 0.27
## Cumulative Var        0.42 0.69
## Proportion Explained  0.61 0.39
## Cumulative Proportion 0.61 1.00
## 
## Mean item complexity =  1.8
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  10  with the objective function =  2.89 with Chi Square =  86729.72
## df of  the model are 1  and the objective function was  0 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0.01 
## 
## The harmonic n.obs is  29999 with the empirical chi square  1.28  with prob <  0.26 
## The total n.obs was  29999  with Likelihood Chi Square =  8.43  with prob <  0.0037 
## 
## Tucker Lewis Index of factoring reliability =  0.999
## RMSEA index =  0.016  and the 90 % confidence intervals are  0.007 0.026
## BIC =  -1.88
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1  ML2
## Correlation of (regression) scores with factors   0.98 0.80
## Multiple R square of scores with factors          0.96 0.64
## Minimum correlation of possible factor scores     0.92 0.27

Rotasi Oblimin

fa_2_oblimin <- fa(
       r = data_final,
       nfactors = 2,
       rotate = "oblimin",
       fm = "ml"
   )
print(fa_2_oblimin)
## Factor Analysis using method =  ml
## Call: fa(r = data_final, nfactors = 2, rotate = "oblimin", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     ML1   ML2   h2    u2 com
## X2 0.75  0.01 0.58 0.423   1
## X3 0.80  0.01 0.66 0.343   1
## X4 0.00  1.00 1.00 0.005   1
## X5 0.78  0.00 0.61 0.386   1
## X6 0.80 -0.02 0.61 0.388   1
## 
##                        ML1  ML2
## SS loadings           2.46 1.00
## Proportion Var        0.49 0.20
## Cumulative Var        0.49 0.69
## Proportion Explained  0.71 0.29
## Cumulative Proportion 0.71 1.00
## 
##  With factor correlations of 
##      ML1  ML2
## ML1 1.00 0.87
## ML2 0.87 1.00
## 
## Mean item complexity =  1
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  10  with the objective function =  2.89 with Chi Square =  86729.72
## df of  the model are 1  and the objective function was  0 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0.01 
## 
## The harmonic n.obs is  29999 with the empirical chi square  1.28  with prob <  0.26 
## The total n.obs was  29999  with Likelihood Chi Square =  8.43  with prob <  0.0037 
## 
## Tucker Lewis Index of factoring reliability =  0.999
## RMSEA index =  0.016  and the 90 % confidence intervals are  0.007 0.026
## BIC =  -1.88
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1  ML2
## Correlation of (regression) scores with factors   0.73 0.80
## Multiple R square of scores with factors          0.53 0.63
## Minimum correlation of possible factor scores     0.06 0.27

Memberi Tanda Loading ≥ 0.4

threshold <- 0.4
annotate_loading <- function(fa_model){
       L <- as.matrix(fa_model$loadings)
       L_annot <- matrix(sprintf("%.3f", L),
                         nrow=nrow(L),
                         dimnames=dimnames(L))
       L_annot[abs(L) >= threshold] <- 
             paste0(L_annot[abs(L) >= threshold], " *")
       return(as.data.frame(L_annot))
   }
cat("\nVarimax Loadings (>=0.4 diberi *):\n")
## 
## Varimax Loadings (>=0.4 diberi *):
print(annotate_loading(fa_2_varimax))
##        ML1     ML2
## X2 0.530 * 0.544 *
## X3 0.567 * 0.579 *
## X4 0.955 *   0.289
## X5 0.545 * 0.563 *
## X6 0.536 * 0.570 *
cat("\nOblimin Loadings (>=0.4 diberi *):\n")
## 
## Oblimin Loadings (>=0.4 diberi *):
print(annotate_loading(fa_2_oblimin))
##        ML1     ML2
## X2 0.753 *   0.008
## X3 0.799 *   0.013
## X4   0.000 0.997 *
## X5 0.781 *   0.003
## X6 0.801 *  -0.021

Pada rotasi varimax, sebagian besar variabel (X2, X3, X5, X6) memiliki loading ≥ 0,4 pada kedua faktor, sehingga masih terjadi cross-loading dan pemisahan faktor belum terlalu jelas, meskipun X4 dominan pada ML1. Sebaliknya, pada rotasi oblimin, struktur faktor menjadi jauh lebih tegas: X2, X3, X5, dan X6 loading tinggi hanya pada ML1, sedangkan X4 sangat dominan pada ML2 (0,997), sehingga setiap variabel lebih jelas terkelompok pada satu faktor utama. Ini menunjukkan bahwa rotasi oblimin menghasilkan interpretasi struktur faktor yang lebih bersih dan mudah dipahami.

Menampilkan Model Varimax Final

fa_2_varimax
## Factor Analysis using method =  ml
## Call: fa(r = data_final, nfactors = 2, rotate = "varimax", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     ML1  ML2   h2    u2 com
## X2 0.53 0.54 0.58 0.423 2.0
## X3 0.57 0.58 0.66 0.343 2.0
## X4 0.95 0.29 1.00 0.005 1.2
## X5 0.54 0.56 0.61 0.386 2.0
## X6 0.54 0.57 0.61 0.388 2.0
## 
##                        ML1  ML2
## SS loadings           2.10 1.36
## Proportion Var        0.42 0.27
## Cumulative Var        0.42 0.69
## Proportion Explained  0.61 0.39
## Cumulative Proportion 0.61 1.00
## 
## Mean item complexity =  1.8
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  10  with the objective function =  2.89 with Chi Square =  86729.72
## df of  the model are 1  and the objective function was  0 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0.01 
## 
## The harmonic n.obs is  29999 with the empirical chi square  1.28  with prob <  0.26 
## The total n.obs was  29999  with Likelihood Chi Square =  8.43  with prob <  0.0037 
## 
## Tucker Lewis Index of factoring reliability =  0.999
## RMSEA index =  0.016  and the 90 % confidence intervals are  0.007 0.026
## BIC =  -1.88
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1  ML2
## Correlation of (regression) scores with factors   0.98 0.80
## Multiple R square of scores with factors          0.96 0.64
## Minimum correlation of possible factor scores     0.92 0.27

Model FA dua faktor dengan rotasi varimax menunjukkan bahwa kedua faktor bersama-sama menjelaskan 69% variasi data (ML1 = 42%, ML2 = 27%), dengan X4 sangat dominan pada ML1, sementara X2, X3, X5, dan X6 memiliki loading cukup tinggi pada kedua faktor (cross-loading). Meskipun masih ada kompleksitas item yang cukup tinggi (mean complexity 1.8), model memiliki goodness-of-fit yang sangat baik (RMSEA 0,016, TLI 0,999), sehingga secara statistik layak digunakan. Model varimax dipilih karena rotasi orthogonal ini menghasilkan distribusi variasi yang lebih seimbang antar faktor dan tetap mempertahankan interpretasi yang stabil dengan kecocokan model yang sangat baik.

Visualisasi

Faktor Analisis

windows(width = 12, height = 8)
fa.diagram(fa_2_varimax, 
                         cut = 0.4,
                         simple = TRUE,
                         cex = 1.2)

Visualisasi menunjukkan hasil analisis faktor dengan dua faktor laten (ML1 dan ML2) yang dihubungkan ke variabel teramati (X2–X6) melalui bobot (loading) yaitu, ML1 kuat pada X4 (1,0) dan X3 (0,6), sedangkan ML2 kuat pada X6 (0,6), X5 (0,6), dan X2 (0,5) dan hanya loading ≥0,4 yang ditampilkan sehingga struktur menjadi lebih sederhana. Sehingga X4 dan X3 mengukur konstruk yang sama (ML1), sementara X6, X5, dan X2 mengukur konstruk lain (ML2), menunjukkan dua dimensi mendasar yang menjelaskan korelasi antar variabel.

Plot Skor Faktor

loadings <- as.data.frame(unclass(fa_2_varimax$loadings))
scores_fa_varimax <- factor.scores(data_final, fa_2_varimax)$scores
plot(scores_fa_varimax[,1],
             scores_fa_varimax[,2],
             xlab = "Faktor 1 (Engagement)",
             ylab = "Faktor 2 (Reach)",
             main = "Plot Skor Faktor",
             pch = 19,
             col = rgb(0,0,1,0.3))

Plot ini menampilkan Faktor 1 (Engagement) pada sumbu horizontal dan Faktor 2 (Reach) pada sumbu vertikal dan setiap titik biru mewakili satu observasi dengan transparansi sehingga area yang padat terlihat lebih gelap, membantu melihat tumpang tindih. Posisi titik menunjukkan skor masing‑masing observasi pada dua faktor tersebut, sehingga titik yang berdekatan memiliki pola respons serupa, sementara titik yang jauh atau terpisah menunjukkan outlier atau kelompok berbeda. Sehingga, interpretasinya berguna untuk mengidentifikasi klaster, hubungan antar faktor, dan observasi yang menonjol.