1.Data Load

Analisis ini menggunakan data Air Quality yang berisi 17 variabel dan 6941 baris. Pada setiap variabel memiliki karateristiknya masing-masing yang dapat dilakukan analisis. Pada proyek ini, data akan dilakukan PCA dan FA dengan tipe data variabel harus numerik/angka sehingga tipe data selain itu akan langsung dilakukan penghapusan. Setelah itu, data harus cek terlebih dahulu bagaimana statistik deskriptifnya, tipe data keseluruhan dan missing value, kemudian dilakukan analisis lebih lanjut.

library(psych)
library(FactoMineR)
library(factoextra)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(ggplot2)
library(corrplot)
## corrplot 0.95 loaded
library(reshape2)
data <- read.csv("AirQualityUCI.csv", sep=";", dec=",")
head(data)
##         Date     Time CO.GT. PT08.S1.CO. NMHC.GT. C6H6.GT. PT08.S2.NMHC.
## 1 10/03/2004 18.00.00    2.6        1360      150     11.9          1046
## 2 10/03/2004 19.00.00    2.0        1292      112      9.4           955
## 3 10/03/2004 20.00.00    2.2        1402       88      9.0           939
## 4 10/03/2004 21.00.00    2.2        1376       80      9.2           948
## 5 10/03/2004 22.00.00    1.6        1272       51      6.5           836
## 6 10/03/2004 23.00.00    1.2        1197       38      4.7           750
##   NOx.GT. PT08.S3.NOx. NO2.GT. PT08.S4.NO2. PT08.S5.O3.    T   RH     AH  X X.1
## 1     166         1056     113         1692        1268 13.6 48.9 0.7578 NA  NA
## 2     103         1174      92         1559         972 13.3 47.7 0.7255 NA  NA
## 3     131         1140     114         1555        1074 11.9 54.0 0.7502 NA  NA
## 4     172         1092     122         1584        1203 11.0 60.0 0.7867 NA  NA
## 5     131         1205     116         1490        1110 11.2 59.6 0.7888 NA  NA
## 6      89         1337      96         1393         949 11.2 59.2 0.7848 NA  NA
summary(data)
##      Date               Time               CO.GT.         PT08.S1.CO.  
##  Length:9471        Length:9471        Min.   :-200.00   Min.   :-200  
##  Class :character   Class :character   1st Qu.:   0.60   1st Qu.: 921  
##  Mode  :character   Mode  :character   Median :   1.50   Median :1053  
##                                        Mean   : -34.21   Mean   :1049  
##                                        3rd Qu.:   2.60   3rd Qu.:1221  
##                                        Max.   :  11.90   Max.   :2040  
##                                        NA's   :114       NA's   :114   
##     NMHC.GT.         C6H6.GT.        PT08.S2.NMHC.       NOx.GT.      
##  Min.   :-200.0   Min.   :-200.000   Min.   :-200.0   Min.   :-200.0  
##  1st Qu.:-200.0   1st Qu.:   4.000   1st Qu.: 711.0   1st Qu.:  50.0  
##  Median :-200.0   Median :   7.900   Median : 895.0   Median : 141.0  
##  Mean   :-159.1   Mean   :   1.866   Mean   : 894.6   Mean   : 168.6  
##  3rd Qu.:-200.0   3rd Qu.:  13.600   3rd Qu.:1105.0   3rd Qu.: 284.0  
##  Max.   :1189.0   Max.   :  63.700   Max.   :2214.0   Max.   :1479.0  
##  NA's   :114      NA's   :114        NA's   :114      NA's   :114     
##   PT08.S3.NOx.     NO2.GT.         PT08.S4.NO2.   PT08.S5.O3.    
##  Min.   :-200   Min.   :-200.00   Min.   :-200   Min.   :-200.0  
##  1st Qu.: 637   1st Qu.:  53.00   1st Qu.:1185   1st Qu.: 700.0  
##  Median : 794   Median :  96.00   Median :1446   Median : 942.0  
##  Mean   : 795   Mean   :  58.15   Mean   :1391   Mean   : 975.1  
##  3rd Qu.: 960   3rd Qu.: 133.00   3rd Qu.:1662   3rd Qu.:1255.0  
##  Max.   :2683   Max.   : 340.00   Max.   :2775   Max.   :2523.0  
##  NA's   :114    NA's   :114       NA's   :114    NA's   :114     
##        T                  RH                AH               X          
##  Min.   :-200.000   Min.   :-200.00   Min.   :-200.0000   Mode:logical  
##  1st Qu.:  10.900   1st Qu.:  34.10   1st Qu.:   0.6923   NA's:9471     
##  Median :  17.200   Median :  48.60   Median :   0.9768                 
##  Mean   :   9.778   Mean   :  39.49   Mean   :  -6.8376                 
##  3rd Qu.:  24.100   3rd Qu.:  61.90   3rd Qu.:   1.2962                 
##  Max.   :  44.600   Max.   :  88.70   Max.   :   2.2310                 
##  NA's   :114        NA's   :114       NA's   :114                       
##    X.1         
##  Mode:logical  
##  NA's:9471     
##                
##                
##                
##                
## 
sapply(data, class)
##          Date          Time        CO.GT.   PT08.S1.CO.      NMHC.GT. 
##   "character"   "character"     "numeric"     "integer"     "integer" 
##      C6H6.GT. PT08.S2.NMHC.       NOx.GT.  PT08.S3.NOx.       NO2.GT. 
##     "numeric"     "integer"     "integer"     "integer"     "integer" 
##  PT08.S4.NO2.   PT08.S5.O3.             T            RH            AH 
##     "integer"     "integer"     "numeric"     "numeric"     "numeric" 
##             X           X.1 
##     "logical"     "logical"
colSums(is.na(data))
##          Date          Time        CO.GT.   PT08.S1.CO.      NMHC.GT. 
##             0             0           114           114           114 
##      C6H6.GT. PT08.S2.NMHC.       NOx.GT.  PT08.S3.NOx.       NO2.GT. 
##           114           114           114           114           114 
##  PT08.S4.NO2.   PT08.S5.O3.             T            RH            AH 
##           114           114           114           114           114 
##             X           X.1 
##          9471          9471
missing_count <- sapply(data, function(x) sum(x == -200, na.rm = TRUE))
missing_count
##          Date          Time        CO.GT.   PT08.S1.CO.      NMHC.GT. 
##             0             0          1683           366          8443 
##      C6H6.GT. PT08.S2.NMHC.       NOx.GT.  PT08.S3.NOx.       NO2.GT. 
##           366           366          1639           366          1642 
##  PT08.S4.NO2.   PT08.S5.O3.             T            RH            AH 
##           366           366           366           366           366 
##             X           X.1 
##             0             0

2. Penanganan Missing Value

Setelah dilakukan pengecekan, terlihat ada missing value pada setiap variabel. Terdapat dua missing value yaitu nilai Null dan -200. Nilai -200 merupakan data gagal terukur oleh sensor, jadi dapat dilakukan konversi menjadi NA yang disebut sebagai missing value sedangkan nilai yang Null ini memiliki jumlah yang sama sebesar 114 sehingga dapat dilakukan penghapusan data. Setelah dilakukan penghapusan, statistik deskriptifnya akan normal dan sesuai.

data_numeric <- data[, !(names(data) %in% c("Date","Time","X","X.1", "NMHC.GT."))]
data_numeric[data_numeric == -200] <- NA
data_now <- na.omit(data_numeric)
colnames(data_now)
##  [1] "CO.GT."        "PT08.S1.CO."   "C6H6.GT."      "PT08.S2.NMHC."
##  [5] "NOx.GT."       "PT08.S3.NOx."  "NO2.GT."       "PT08.S4.NO2." 
##  [9] "PT08.S5.O3."   "T"             "RH"            "AH"
summary(data_now)
##      CO.GT.        PT08.S1.CO.      C6H6.GT.     PT08.S2.NMHC.   
##  Min.   : 0.100   Min.   : 647   Min.   : 0.20   Min.   : 390.0  
##  1st Qu.: 1.100   1st Qu.: 956   1st Qu.: 4.90   1st Qu.: 760.0  
##  Median : 1.900   Median :1085   Median : 8.80   Median : 931.0  
##  Mean   : 2.182   Mean   :1120   Mean   :10.55   Mean   : 958.5  
##  3rd Qu.: 2.900   3rd Qu.:1254   3rd Qu.:14.60   3rd Qu.:1135.0  
##  Max.   :11.900   Max.   :2040   Max.   :63.70   Max.   :2214.0  
##     NOx.GT.        PT08.S3.NOx.       NO2.GT.       PT08.S4.NO2. 
##  Min.   :   2.0   Min.   : 322.0   Min.   :  2.0   Min.   : 551  
##  1st Qu.: 103.0   1st Qu.: 642.0   1st Qu.: 79.0   1st Qu.:1207  
##  Median : 186.0   Median : 786.0   Median :110.0   Median :1457  
##  Mean   : 250.7   Mean   : 816.9   Mean   :113.9   Mean   :1453  
##  3rd Qu.: 335.0   3rd Qu.: 947.0   3rd Qu.:142.0   3rd Qu.:1683  
##  Max.   :1479.0   Max.   :2683.0   Max.   :333.0   Max.   :2775  
##   PT08.S5.O3.         T               RH              AH        
##  Min.   : 221   Min.   :-1.90   Min.   : 9.20   Min.   :0.1847  
##  1st Qu.: 760   1st Qu.:11.20   1st Qu.:35.30   1st Qu.:0.6941  
##  Median :1006   Median :16.80   Median :49.20   Median :0.9539  
##  Mean   :1058   Mean   :17.76   Mean   :48.88   Mean   :0.9856  
##  3rd Qu.:1322   3rd Qu.:23.70   3rd Qu.:62.20   3rd Qu.:1.2516  
##  Max.   :2523   Max.   :44.60   Max.   :88.70   Max.   :2.1806
par(mfrow=c(2,3))
for(i in 1:ncol(data_now)){
  hist(data_now[,i],
       main=colnames(data_now)[i],
       col="lightblue")
}

3. Pengecekan Outlier

Setelah penanganan missing value, dilakukan deteksi outlier menggunakan metode Z-score dengan kriteria |Z| > 3. Pemeriksaan ini bertujuan untuk mengidentifikasi nilai ekstrem yang berpotensi memengaruhi matriks korelasi dalam analisis PCA dan faktor analisis. Hasil menunjukkan bahwa proporsi outlier pada masing-masing variabel relatif kecil (< 5%), sehingga data masih dianggap representatif dan tidak dilakukan penghapusan.

z_scores      <- scale(data_now)
outlier_count <- colSums(abs(z_scores) > 3)
print(data.frame(Variable = names(outlier_count), Outliers_Z3 = outlier_count))
##                    Variable Outliers_Z3
## CO.GT.               CO.GT.          87
## PT08.S1.CO.     PT08.S1.CO.          41
## C6H6.GT.           C6H6.GT.          85
## PT08.S2.NMHC. PT08.S2.NMHC.          27
## NOx.GT.             NOx.GT.         118
## PT08.S3.NOx.   PT08.S3.NOx.          87
## NO2.GT.             NO2.GT.          36
## PT08.S4.NO2.   PT08.S4.NO2.          27
## PT08.S5.O3.     PT08.S5.O3.          27
## T                         T           2
## RH                       RH           0
## AH                       AH           0
n_rows_plot <- 2
n_cols_plot <- 5

colors <- rainbow(ncol(data_now))

par(mfrow = c(n_rows_plot, n_cols_plot), mar = c(4, 4, 3, 1))

for (i in seq_along(colnames(data_now))) {
  col <- colnames(data_now)[i]
  boxplot(data_now[[col]],
          main       = col,
          col        = colors[i],
          border     = "gray30",
          horizontal = TRUE)
}

4.Correlation Matriks

cor_matrix <- cor(data_now)
round(cor_matrix, 3)
##               CO.GT. PT08.S1.CO. C6H6.GT. PT08.S2.NMHC. NOx.GT. PT08.S3.NOx.
## CO.GT.         1.000       0.877    0.930         0.914   0.786       -0.701
## PT08.S1.CO.    0.877       1.000    0.877         0.886   0.708       -0.763
## C6H6.GT.       0.930       0.877    1.000         0.983   0.718       -0.726
## PT08.S2.NMHC.  0.914       0.886    0.983         1.000   0.705       -0.782
## NOx.GT.        0.786       0.708    0.718         0.705   1.000       -0.662
## PT08.S3.NOx.  -0.701      -0.763   -0.726        -0.782  -0.662        1.000
## NO2.GT.        0.674       0.628    0.603         0.633   0.757       -0.641
## PT08.S4.NO2.   0.631       0.676    0.762         0.774   0.234       -0.511
## PT08.S5.O3.    0.853       0.897    0.861         0.877   0.789       -0.793
## T              0.018       0.028    0.189         0.228  -0.276       -0.099
## RH             0.065       0.169   -0.022        -0.046   0.232       -0.116
## AH             0.059       0.150    0.187         0.206  -0.144       -0.223
##               NO2.GT. PT08.S4.NO2. PT08.S5.O3.      T     RH     AH
## CO.GT.          0.674        0.631       0.853  0.018  0.065  0.059
## PT08.S1.CO.     0.628        0.676       0.897  0.028  0.169  0.150
## C6H6.GT.        0.603        0.762       0.861  0.189 -0.022  0.187
## PT08.S2.NMHC.   0.633        0.774       0.877  0.228 -0.046  0.206
## NOx.GT.         0.757        0.234       0.789 -0.276  0.232 -0.144
## PT08.S3.NOx.   -0.641       -0.511      -0.793 -0.099 -0.116 -0.223
## NO2.GT.         1.000        0.143       0.703 -0.214 -0.075 -0.350
## PT08.S4.NO2.    0.143        1.000       0.574  0.567 -0.009  0.646
## PT08.S5.O3.     0.703        0.574       1.000 -0.046  0.165  0.076
## T              -0.214        0.567      -0.046  1.000 -0.564  0.661
## RH             -0.075       -0.009       0.165 -0.564  1.000  0.180
## AH             -0.350        0.646       0.076  0.661  0.180  1.000
cor_matrix <- cor(data_now)
corrplot(cor_matrix,
         method = "color",
         type = "upper",
         addCoef.col = "black",
         tl.col = "black",
         tl.srt = 45,
         col = colorRampPalette(c("blue", "white", "red"))(200))

5. UJI ASUMSI Kaiser-Meyer-Olkin (KMO) & Bartlett

Pada uji asumsi terdapat dua uji yaitu KMO dan Bartlett. Kedua uji ini memiliki fungsi untuk menilai kelayakan data sebelum dilakukan analisis faktor. Uji KMO digunakan untuk mengukur kecukupan sampel dan melihat apakah pola korelasi antar variabel cukup kompak untuk membentuk faktor. Nilai KMO berada pada rentang 0 - 1, jika nilainya > 0,5 akan dikatakan layak sedangkan < 0,5 tidak layak/gunakan. Bartlett’s Test of Sphericity digunakan untuk menguji apakah matriks korelasi berbeda secara signifikan dari matriks identitas. Jika nilai signifikansi (p-value) < 0,05, maka terdapat korelasi yang signifikan antar variabel, sehingga analisis faktor dapat dilanjutkan.

kmo_result <- KMO(data_now)

interpret_kmo <- function(val) {
  if      (val >= 0.90) "Marvelous"
  else if (val >= 0.80) "Meritorious"
  else if (val >= 0.70) "Middling"
  else if (val >= 0.60) "Mediocre"
  else                  "Unacceptable"
}

cat("Overall KMO    :", round(kmo_result$MSA, 4), "\n")
## Overall KMO    : 0.817
cat("Interpretation :", interpret_kmo(kmo_result$MSA), "\n\n")
## Interpretation : Meritorious
# Hitung KMO
kmo_result <- KMO(data_now)

# Lihat keseluruhan
kmo_result$MSA
## [1] 0.8170217
# Lihat per variabel
kmo_result$MSAi
##        CO.GT.   PT08.S1.CO.      C6H6.GT. PT08.S2.NMHC.       NOx.GT. 
##     0.9509872     0.9352423     0.8613277     0.8342850     0.8420433 
##  PT08.S3.NOx.       NO2.GT.  PT08.S4.NO2.   PT08.S5.O3.             T 
##     0.8508353     0.9194885     0.7758389     0.9489606     0.4701238 
##            RH            AH 
##     0.2530623     0.4640117
cor_matrix <- cor(data_now)
cortest.bartlett(cor_matrix, n = nrow(data_now))
## $chisq
## [1] 134460.3
## 
## $p.value
## [1] 0
## 
## $df
## [1] 66

Pada hasil diatas terlihat bahwa masih ada nilai KMO < 0,5 yaitu T, RH, dan AH. Ketiga variabel ini akan dilakukan penghapusan secara satu-persatu. Dari ketiga nilai, kita ambil nilai paling kecil yaitu RH dengan nilai 0.25, kemudian jika sudah dihapus nilai KMO pada masing-masing variabel akn berubah/semakin naik.

data_step1 <- subset(data_now, select = -RH)
data_clean <- na.omit(data_step1)
head(data_clean, 10)
##    CO.GT. PT08.S1.CO. C6H6.GT. PT08.S2.NMHC. NOx.GT. PT08.S3.NOx. NO2.GT.
## 1     2.6        1360     11.9          1046     166         1056     113
## 2     2.0        1292      9.4           955     103         1174      92
## 3     2.2        1402      9.0           939     131         1140     114
## 4     2.2        1376      9.2           948     172         1092     122
## 5     1.6        1272      6.5           836     131         1205     116
## 6     1.2        1197      4.7           750      89         1337      96
## 7     1.2        1185      3.6           690      62         1462      77
## 8     1.0        1136      3.3           672      62         1453      76
## 9     0.9        1094      2.3           609      45         1579      60
## 12    0.7        1066      1.1           512      16         1918      28
##    PT08.S4.NO2. PT08.S5.O3.    T     AH
## 1          1692        1268 13.6 0.7578
## 2          1559         972 13.3 0.7255
## 3          1555        1074 11.9 0.7502
## 4          1584        1203 11.0 0.7867
## 5          1490        1110 11.2 0.7888
## 6          1393         949 11.2 0.7848
## 7          1333         733 11.3 0.7603
## 8          1333         730 10.7 0.7702
## 9          1276         620 10.7 0.7648
## 12         1182         422 11.0 0.7366
colSums(is.na(data_clean))
##        CO.GT.   PT08.S1.CO.      C6H6.GT. PT08.S2.NMHC.       NOx.GT. 
##             0             0             0             0             0 
##  PT08.S3.NOx.       NO2.GT.  PT08.S4.NO2.   PT08.S5.O3.             T 
##             0             0             0             0             0 
##            AH 
##             0
cor_matrix1 <- cor(data_step1)
KMO(cor_matrix1)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor_matrix1)
## Overall MSA =  0.84
## MSA for each item = 
##        CO.GT.   PT08.S1.CO.      C6H6.GT. PT08.S2.NMHC.       NOx.GT. 
##          0.95          0.93          0.87          0.81          0.82 
##  PT08.S3.NOx.       NO2.GT.  PT08.S4.NO2.   PT08.S5.O3.             T 
##          0.82          0.89          0.74          0.94          0.64 
##            AH 
##          0.45
cortest.bartlett(cor_matrix)
## Warning in cortest.bartlett(cor_matrix): n not specified, 100 used
## $chisq
## [1] 1825.72
## 
## $p.value
## [1] 0
## 
## $df
## [1] 66

Setelah dilakukan penghapusan terhadap variabel RH, terlihat ada 1 variabel yang masih memiliki nilai KMO < 0.5 yaitu variabel AH, maka dilakukan penghapusan kembali.

data_step2 <- subset(data_clean, select = -AH)
head(data_step2, 10)
##    CO.GT. PT08.S1.CO. C6H6.GT. PT08.S2.NMHC. NOx.GT. PT08.S3.NOx. NO2.GT.
## 1     2.6        1360     11.9          1046     166         1056     113
## 2     2.0        1292      9.4           955     103         1174      92
## 3     2.2        1402      9.0           939     131         1140     114
## 4     2.2        1376      9.2           948     172         1092     122
## 5     1.6        1272      6.5           836     131         1205     116
## 6     1.2        1197      4.7           750      89         1337      96
## 7     1.2        1185      3.6           690      62         1462      77
## 8     1.0        1136      3.3           672      62         1453      76
## 9     0.9        1094      2.3           609      45         1579      60
## 12    0.7        1066      1.1           512      16         1918      28
##    PT08.S4.NO2. PT08.S5.O3.    T
## 1          1692        1268 13.6
## 2          1559         972 13.3
## 3          1555        1074 11.9
## 4          1584        1203 11.0
## 5          1490        1110 11.2
## 6          1393         949 11.2
## 7          1333         733 11.3
## 8          1333         730 10.7
## 9          1276         620 10.7
## 12         1182         422 11.0

Step 2 pada penghapusan variabel AH telah berhasil sehingga nilai KMO sudah > 0.5 dengan 10 variabel dan hasil Bartlett’s Test menunjukkan nilai χ² sebesar 1465.287 dengan p-value < 0.05 (1.116829e-277), sehingga H0 ditolak. Hal ini menunjukkan bahwa terdapat korelasi signifikan antar variabel dan data layak untuk dilakukan analisis faktor.

cor_matrix2 <- cor(data_step2)
KMO(cor_matrix2)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor_matrix2)
## Overall MSA =  0.88
## MSA for each item = 
##        CO.GT.   PT08.S1.CO.      C6H6.GT. PT08.S2.NMHC.       NOx.GT. 
##          0.95          0.93          0.85          0.84          0.92 
##  PT08.S3.NOx.       NO2.GT.  PT08.S4.NO2.   PT08.S5.O3.             T 
##          0.91          0.87          0.79          0.94          0.56
cortest.bartlett(cor_matrix2)
## Warning in cortest.bartlett(cor_matrix2): n not specified, 100 used
## $chisq
## [1] 1465.287
## 
## $p.value
## [1] 1.116829e-277
## 
## $df
## [1] 45

6. Data Scaled

Setelah uji asumsi (KMO dan Bartlett) menunjukkan bahwa data layak untuk Analisis Faktor, langkah selanjutnya adalah melakukan standardisasi terhadap data. Scaling dilakukan karena variabel memiliki satuan dan rentang nilai yang berbeda. Jika tidak distandarisasi, variabel dengan skala besar akan lebih dominan dalam pembentukan faktor.Setelah distandarisasi, data siap digunakan untuk proses ekstraksi faktor.

data_scaled <- scale(data_step2)
head(data_scaled, 10)
##         CO.GT. PT08.S1.CO.   C6H6.GT. PT08.S2.NMHC.    NOx.GT. PT08.S3.NOx.
## 1   0.28972086   1.0976209  0.1802373    0.33120749 -0.4058837     0.949223
## 2  -0.12661104   0.7867406 -0.1546487   -0.01341768 -0.7078806     1.417668
## 3   0.01216626   1.2896351 -0.2082305   -0.07401112 -0.5736598     1.282692
## 4   0.01216626   1.1707692 -0.1814396   -0.03992731 -0.3771221     1.092138
## 5  -0.40416564   0.6953053 -0.5431165   -0.46408136 -0.5736598     1.540734
## 6  -0.68172024   0.3524227 -0.7842345   -0.78977109 -0.7749911     2.064757
## 7  -0.68172024   0.2975614 -0.9315843   -1.01699647 -0.9044183     2.560992
## 8  -0.82049754   0.0735448 -0.9717706   -1.08516409 -0.9044183     2.525263
## 9  -0.88988619  -0.1184695 -1.1057251   -1.32375075 -0.9859096     3.025467
## 12 -1.02866349  -0.2464790 -1.2664704   -1.69109846 -1.1249241     4.371254
##        NO2.GT. PT08.S4.NO2. PT08.S5.O3.          T
## 1  -0.01841140    0.6774697  0.51719147 -0.4697983
## 2  -0.46074932    0.3010207 -0.21095798 -0.5037161
## 3   0.00265231    0.2896990  0.03995838 -0.6619993
## 4   0.17116199    0.3717818  0.35729378 -0.7637527
## 5   0.04477973    0.1057202  0.12851710 -0.7411409
## 6  -0.37649447   -0.1688328 -0.26753716 -0.7411409
## 7  -0.77670497   -0.3386594 -0.79888946 -0.7298349
## 8  -0.79776868   -0.3386594 -0.80626935 -0.7976706
## 9  -1.13478804   -0.4999947 -1.07686543 -0.7976706
## 12 -1.80882677   -0.7660563 -1.56393837 -0.7637527

#7. Pembentukan Principal Component Analysis Berdasarkan kriteria eigenvalue > 1 dan nilai cumulative proportion sebesar 86.28%, dapat disimpulkan bahwa 2 komponen utama sudah cukup untuk mewakili keseluruhan variabel dalam dataset. Besarnya kontribusi PC1 (68.81%) menunjukkan bahwa terdapat satu komponen dominan yang sangat kuat dalam data. Hal ini sejalan dengan hasil korelasi dan analisis faktor sebelumnya, di mana variabel-variabel polutan memiliki hubungan yang sangat tinggi satu sama lain. PC2 memberikan tambahan variasi sebesar 17.47%, yang kemungkinan merepresentasikan pola kedua dalam data, seperti pengaruh suhu atau karakteristik lingkungan.

pca_result <- prcomp(data_scaled, scale. = FALSE)
summary(pca_result)
## Importance of components:
##                           PC1    PC2     PC3     PC4    PC5     PC6     PC7
## Standard deviation     2.6232 1.3216 0.68793 0.57205 0.4615 0.35743 0.30970
## Proportion of Variance 0.6881 0.1747 0.04732 0.03272 0.0213 0.01278 0.00959
## Cumulative Proportion  0.6881 0.8628 0.91009 0.94282 0.9641 0.97689 0.98648
##                            PC8    PC9    PC10
## Standard deviation     0.26587 0.2324 0.10251
## Proportion of Variance 0.00707 0.0054 0.00105
## Cumulative Proportion  0.99355 0.9990 1.00000
eigenvalues <- pca_result$sdev^2
eigenvalues
##  [1] 6.88105245 1.74662722 0.47324440 0.32723929 0.21297519 0.12775749
##  [7] 0.09591667 0.07068543 0.05399304 0.01050882
plot(eigenvalues, type="b",
     xlab="Komponen",
     ylab="Eigenvalue",
     main="Scree Plot")
abline(h=1, col="red", lty=2)

round(pca_result$rotation, 3)
##                  PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8    PC9
## CO.GT.        -0.360 -0.026 -0.223 -0.281  0.053  0.367 -0.310  0.470 -0.533
## PT08.S1.CO.   -0.357  0.010 -0.215  0.225 -0.297 -0.331 -0.688 -0.041  0.320
## C6H6.GT.      -0.365  0.116 -0.166 -0.215  0.079  0.228  0.241  0.107  0.478
## PT08.S2.NMHC. -0.369  0.131 -0.034 -0.115  0.004  0.188  0.265  0.044  0.406
## NOx.GT.       -0.309 -0.338  0.025 -0.230  0.701 -0.255 -0.111 -0.404 -0.058
## PT08.S3.NOx.   0.320  0.018 -0.480 -0.699 -0.199 -0.357  0.039 -0.024  0.051
## NO2.GT.       -0.279 -0.350  0.535 -0.369 -0.553  0.036  0.043 -0.235 -0.093
## PT08.S4.NO2.  -0.261  0.508 -0.277  0.089 -0.166  0.075  0.152 -0.620 -0.384
## PT08.S5.O3.   -0.358 -0.087 -0.070  0.217 -0.088 -0.621  0.480  0.360 -0.241
## T             -0.032  0.683  0.525 -0.278  0.174 -0.290 -0.180  0.178 -0.007
##                 PC10
## CO.GT.        -0.052
## PT08.S1.CO.   -0.004
## C6H6.GT.       0.651
## PT08.S2.NMHC. -0.748
## NOx.GT.       -0.028
## PT08.S3.NOx.  -0.076
## NO2.GT.        0.066
## PT08.S4.NO2.   0.043
## PT08.S5.O3.    0.036
## T              0.025
loadings <- pca_result$rotation
round(loadings, 3)
##                  PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8    PC9
## CO.GT.        -0.360 -0.026 -0.223 -0.281  0.053  0.367 -0.310  0.470 -0.533
## PT08.S1.CO.   -0.357  0.010 -0.215  0.225 -0.297 -0.331 -0.688 -0.041  0.320
## C6H6.GT.      -0.365  0.116 -0.166 -0.215  0.079  0.228  0.241  0.107  0.478
## PT08.S2.NMHC. -0.369  0.131 -0.034 -0.115  0.004  0.188  0.265  0.044  0.406
## NOx.GT.       -0.309 -0.338  0.025 -0.230  0.701 -0.255 -0.111 -0.404 -0.058
## PT08.S3.NOx.   0.320  0.018 -0.480 -0.699 -0.199 -0.357  0.039 -0.024  0.051
## NO2.GT.       -0.279 -0.350  0.535 -0.369 -0.553  0.036  0.043 -0.235 -0.093
## PT08.S4.NO2.  -0.261  0.508 -0.277  0.089 -0.166  0.075  0.152 -0.620 -0.384
## PT08.S5.O3.   -0.358 -0.087 -0.070  0.217 -0.088 -0.621  0.480  0.360 -0.241
## T             -0.032  0.683  0.525 -0.278  0.174 -0.290 -0.180  0.178 -0.007
##                 PC10
## CO.GT.        -0.052
## PT08.S1.CO.   -0.004
## C6H6.GT.       0.651
## PT08.S2.NMHC. -0.748
## NOx.GT.       -0.028
## PT08.S3.NOx.  -0.076
## NO2.GT.        0.066
## PT08.S4.NO2.   0.043
## PT08.S5.O3.    0.036
## T              0.025
eig_tabel <- data.frame(
  PC         = paste0("PC", 1:length(eigenvalues)),
  Eigenvalue = round(eigenvalues, 4),
  Pct_Var    = round(eigenvalues / sum(eigenvalues) * 100, 2),
  Cum_Var    = round(cumsum(eigenvalues / sum(eigenvalues) * 100), 2)
)
print(eig_tabel)
##      PC Eigenvalue Pct_Var Cum_Var
## 1   PC1     6.8811   68.81   68.81
## 2   PC2     1.7466   17.47   86.28
## 3   PC3     0.4732    4.73   91.01
## 4   PC4     0.3272    3.27   94.28
## 5   PC5     0.2130    2.13   96.41
## 6   PC6     0.1278    1.28   97.69
## 7   PC7     0.0959    0.96   98.65
## 8   PC8     0.0707    0.71   99.35
## 9   PC9     0.0540    0.54   99.89
## 10 PC10     0.0105    0.11  100.00
n_factors <- sum(eigenvalues > 1)
cat("Jumlah faktor dipertahankan (Kaiser Rule):", n_factors, "\n")
## Jumlah faktor dipertahankan (Kaiser Rule): 2
cat("Kumulatif varians yang dijelaskan        :", eig_tabel$Cum_Var[n_factors], "%\n")
## Kumulatif varians yang dijelaskan        : 86.28 %

8. Factor Analysis

FA - Stage 5a: Unrotated Factor Matrix

FA dijalankan tanpa rotasi terlebih dahulu untuk melihat pola loading awal dan komunalitas setiap variabel sebelum rotasi dilakukan. Komunalitas < 0.50 menandakan variabel kurang berkontribusi pada faktor.

fa_unrotated <- fa(r        = cor_matrix2,
                   nfactors = n_factors,
                   rotate   = "none",
                   fm       = "pa",
                   n.obs    = nrow(data_step2))

cat("=== Factor Loadings UNROTATED (cutoff = 0.40) ===\n")
## === Factor Loadings UNROTATED (cutoff = 0.40) ===
print(fa_unrotated$loadings, cutoff = 0.40, digits = 3)
## 
## Loadings:
##               PA1    PA2   
## CO.GT.         0.937       
## PT08.S1.CO.    0.927       
## C6H6.GT.       0.960       
## PT08.S2.NMHC.  0.977       
## NOx.GT.        0.798 -0.454
## PT08.S3.NOx.  -0.805       
## NO2.GT.        0.704 -0.427
## PT08.S4.NO2.   0.688  0.701
## PT08.S5.O3.    0.934       
## T                     0.733
## 
##                  PA1   PA2
## SS loadings    6.745 1.487
## Proportion Var 0.675 0.149
## Cumulative Var 0.675 0.823
cat("=== Komunalitas Unrotated ===\n")
## === Komunalitas Unrotated ===
kom_unrot <- round(fa_unrotated$communality, 3)
print(kom_unrot)
##        CO.GT.   PT08.S1.CO.      C6H6.GT. PT08.S2.NMHC.       NOx.GT. 
##         0.881         0.859         0.944         0.984         0.842 
##  PT08.S3.NOx.       NO2.GT.  PT08.S4.NO2.   PT08.S5.O3.             T 
##         0.649         0.677         0.964         0.888         0.544
var_rendah <- names(kom_unrot[kom_unrot < 0.50])
if (length(var_rendah) > 0) {
  cat("Variabel dengan komunalitas < 0.50:", var_rendah, "\n")
} else {
  cat("Semua komunalitas >= 0.50 -> lanjut ke rotasi.\n")
}
## Semua komunalitas >= 0.50 -> lanjut ke rotasi.

FA - Stage 5b: Rotated Factor Matrix (Varimax)

Rotasi Varimax diterapkan agar setiap variabel memiliki loading tinggi pada satu faktor saja sehingga interpretasi lebih mudah. Hasil ini merupakan hasil utama FA yang digunakan untuk menamai dan menginterpretasikan faktor.

fa_rotated <- fa(r        = cor_matrix2,
                 nfactors = n_factors,
                 rotate   = "varimax",
                 fm       = "pa",
                 n.obs    = nrow(data_step2))

cat("=== Factor Loadings ROTATED - Varimax (cutoff = 0.40) ===\n")
## === Factor Loadings ROTATED - Varimax (cutoff = 0.40) ===
print(fa_rotated$loadings, cutoff = 0.40, digits = 3)
## 
## Loadings:
##               PA1    PA2   
## CO.GT.         0.918       
## PT08.S1.CO.    0.896       
## C6H6.GT.       0.892       
## PT08.S2.NMHC.  0.902  0.412
## NOx.GT.        0.887       
## PT08.S3.NOx.  -0.788       
## NO2.GT.        0.789       
## PT08.S4.NO2.   0.488  0.852
## PT08.S5.O3.    0.936       
## T                     0.730
## 
##                  PA1   PA2
## SS loadings    6.409 1.823
## Proportion Var 0.641 0.182
## Cumulative Var 0.641 0.823
cat("=== Komunalitas Rotated ===\n")
## === Komunalitas Rotated ===
round(fa_rotated$communality, 3)
##        CO.GT.   PT08.S1.CO.      C6H6.GT. PT08.S2.NMHC.       NOx.GT. 
##         0.881         0.859         0.944         0.984         0.842 
##  PT08.S3.NOx.       NO2.GT.  PT08.S4.NO2.   PT08.S5.O3.             T 
##         0.649         0.677         0.964         0.888         0.544
cat("=== Proporsi Varians Tiap Faktor (Varimax) ===\n")
## === Proporsi Varians Tiap Faktor (Varimax) ===
var_acc  <- fa_rotated$Vaccounted
var_tabel <- data.frame(
  SS_Loadings = round(var_acc["SS loadings",   ], 4),
  Prop_Var    = round(var_acc["Proportion Var", ], 4),
  Cum_Var     = round(var_acc["Cumulative Var", ], 4)
)
print(var_tabel)
##     SS_Loadings Prop_Var Cum_Var
## PA1      6.4090   0.6409  0.6409
## PA2      1.8232   0.1823  0.8232

FA - Stage 5b: Visualisasi (Factor Diagram & Biplot)

Factor diagram menampilkan hubungan antar variabel dan faktor secara visual. Biplot menggambarkan arah dan kekuatan loading tiap variabel pada ruang dua faktor.

fa.diagram(fa_rotated,
           main   = "Factor Diagram - Air Quality (10 Variabel, Varimax)",
           cut    = 0.40,
           digits = 2)

# Biplot FA
loading_df <- data.frame(
  variabel = rownames(fa_rotated$loadings),
  FA1      = as.numeric(fa_rotated$loadings[, 1]),
  FA2      = as.numeric(fa_rotated$loadings[, 2])
)

ggplot(loading_df, aes(x = FA1, y = FA2, label = variabel)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray60") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
  geom_segment(aes(x = 0, y = 0, xend = FA1, yend = FA2),
               arrow    = arrow(length = unit(0.2, "cm")),
               color    = "#2E4057", linewidth = 0.7) +
  geom_text(vjust = -0.5, hjust = 0.5, size = 3.5, color = "#1F6B75") +
  xlim(-1.1, 1.1) + ylim(-1.1, 1.1) +
  labs(
    title = "Biplot Factor Analysis (Varimax) - Faktor 1 vs Faktor 2",
    x     = paste0("Faktor 1 (", round(var_acc["Proportion Var", 1] * 100, 1), "%)"),
    y     = paste0("Faktor 2 (", round(var_acc["Proportion Var", 2] * 100, 1), "%)")
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

## FA - Stage 6: Validasi Split Sample Data dibagi 50:50 secara acak untuk menguji apakah struktur faktor yang terbentuk konsisten pada dua subset independen. Jika pola loading kedua sampel mirip, struktur faktor dinyatakan valid dan stabil.

set.seed(123)
idx     <- sample(1:nrow(data_step2), size = floor(nrow(data_step2) / 2))
data_s1 <- data_step2[ idx, ]
data_s2 <- data_step2[-idx, ]

fa_s1 <- fa(r = cor(data_s1), nfactors = n_factors,
            rotate = "varimax", fm = "pa", n.obs = nrow(data_s1))
fa_s2 <- fa(r = cor(data_s2), nfactors = n_factors,
            rotate = "varimax", fm = "pa", n.obs = nrow(data_s2))

cat("=== Loadings Sample 1 (cutoff 0.40) ===\n")
## === Loadings Sample 1 (cutoff 0.40) ===
print(fa_s1$loadings, cutoff = 0.40, digits = 3)
## 
## Loadings:
##               PA1    PA2   
## CO.GT.         0.916       
## PT08.S1.CO.    0.897       
## C6H6.GT.       0.885       
## PT08.S2.NMHC.  0.897  0.423
## NOx.GT.        0.892       
## PT08.S3.NOx.  -0.786       
## NO2.GT.        0.784       
## PT08.S4.NO2.   0.490  0.853
## PT08.S5.O3.    0.935       
## T                     0.732
## 
##                  PA1   PA2
## SS loadings    6.385 1.849
## Proportion Var 0.639 0.185
## Cumulative Var 0.639 0.823
cat("\n=== Loadings Sample 2 (cutoff 0.40) ===\n")
## 
## === Loadings Sample 2 (cutoff 0.40) ===
print(fa_s2$loadings, cutoff = 0.40, digits = 3)
## 
## Loadings:
##               PA1    PA2   
## CO.GT.         0.921       
## PT08.S1.CO.    0.895       
## C6H6.GT.       0.898       
## PT08.S2.NMHC.  0.907  0.401
## NOx.GT.        0.881       
## PT08.S3.NOx.  -0.790       
## NO2.GT.        0.793       
## PT08.S4.NO2.   0.487  0.851
## PT08.S5.O3.    0.937       
## T                     0.728
## 
##                  PA1   PA2
## SS loadings    6.433 1.799
## Proportion Var 0.643 0.180
## Cumulative Var 0.643 0.823

FA - Stage 7: Skor Faktor

Skor faktor adalah nilai setiap observasi pada masing-masing faktor laten. Skor ini bisa digunakan sebagai variabel baru untuk analisis lanjutan seperti clustering atau regresi.

fa_scores   <- factor.scores(x = data_step2, f = fa_rotated)
skor_faktor <- as.data.frame(fa_scores$scores)
colnames(skor_faktor) <- c("Faktor1_Polutan_Pembakaran", "Faktor2_Suhu_NO2")

cat("=== 6 Baris Pertama Skor Faktor ===\n")
## === 6 Baris Pertama Skor Faktor ===
head(round(skor_faktor, 4))
##   Faktor1_Polutan_Pembakaran Faktor2_Suhu_NO2
## 1                     0.3373           0.4970
## 2                    -0.0635           0.4004
## 3                     0.1059           0.1115
## 4                     0.1798           0.1227
## 5                    -0.2090          -0.0250
## 6                    -0.5347          -0.1714
cat("\n=== Statistika Deskriptif Skor Faktor ===\n")
## 
## === Statistika Deskriptif Skor Faktor ===
summary(skor_faktor)
##  Faktor1_Polutan_Pembakaran Faktor2_Suhu_NO2 
##  Min.   :-1.9834            Min.   :-2.1730  
##  1st Qu.:-0.7576            1st Qu.:-0.8490  
##  Median :-0.1335            Median : 0.1121  
##  Mean   : 0.0000            Mean   : 0.0000  
##  3rd Qu.: 0.6106            3rd Qu.: 0.7952  
##  Max.   : 4.3345            Max.   : 3.5528