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
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")
}
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)
}
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))
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
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 %
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.
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
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
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