Implementasi Principal Component Analysis (PCA) dan Factor Analysis (FA) pada Data Kualitas Udara di Kota Italia
Refochella Florentine Gwacevo Pertiwi (24031554045) & Gita Nurani (24031554122)
2026-03-03
library(readxl)
library(psych)
library(GPArotation)##
## Attaching package: 'GPArotation'## The following objects are masked from 'package:psych':
##
## equamax, variminlibrary(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(knitr)
library(corrplot)## corrplot 0.95 loaded##Import Dataset Dataset AirQualityUCI.xlsx berisi data kualitas udara di sebuah kota di Italia dengan berbagai parameter polutan dan sensor pendukung. Data diimpor menggunakan read_excel() untuk selanjutnya dilakukan proses pembersihan dan analisis.
data <- read_excel("AirQualityUCI.xlsx")
str(data)## tibble [9,357 × 15] (S3: tbl_df/tbl/data.frame)
## $ Date : POSIXct[1:9357], format: "2004-03-10" "2004-03-10" ...
## $ Time : POSIXct[1:9357], format: "1899-12-31 18:00:00" "1899-12-31 19:00:00" ...
## $ CO(GT) : num [1:9357] 2.6 2 2.2 2.2 1.6 1.2 1.2 1 0.9 0.6 ...
## $ PT08.S1(CO) : num [1:9357] 1360 1292 1402 1376 1272 ...
## $ NMHC(GT) : num [1:9357] 150 112 88 80 51 38 31 31 24 19 ...
## $ C6H6(GT) : num [1:9357] 11.88 9.4 9 9.23 6.52 ...
## $ PT08.S2(NMHC): num [1:9357] 1046 955 939 948 836 ...
## $ NOx(GT) : num [1:9357] 166 103 131 172 131 89 62 62 45 -200 ...
## $ PT08.S3(NOx) : num [1:9357] 1056 1174 1140 1092 1205 ...
## $ NO2(GT) : num [1:9357] 113 92 114 122 116 96 77 76 60 -200 ...
## $ PT08.S4(NO2) : num [1:9357] 1692 1559 1554 1584 1490 ...
## $ PT08.S5(O3) : num [1:9357] 1268 972 1074 1203 1110 ...
## $ T : num [1:9357] 13.6 13.3 11.9 11 11.2 ...
## $ RH : num [1:9357] 48.9 47.7 54 60 59.6 ...
## $ AH : num [1:9357] 0.758 0.725 0.75 0.787 0.789 ...summary(data)## Date Time CO(GT)
## Min. :2004-03-10 00:00:00 Min. :1899-12-31 00:00:00 Min. :-200.00
## 1st Qu.:2004-06-16 00:00:00 1st Qu.:1899-12-31 05:00:00 1st Qu.: 0.60
## Median :2004-09-21 00:00:00 Median :1899-12-31 11:00:00 Median : 1.50
## Mean :2004-09-21 04:30:05 Mean :1899-12-31 11:29:54 Mean : -34.21
## 3rd Qu.:2004-12-28 00:00:00 3rd Qu.:1899-12-31 18:00:00 3rd Qu.: 2.60
## Max. :2005-04-04 00:00:00 Max. :1899-12-31 23:00:00 Max. : 11.90
## PT08.S1(CO) NMHC(GT) C6H6(GT) PT08.S2(NMHC)
## Min. :-200 Min. :-200.0 Min. :-200.000 Min. :-200.0
## 1st Qu.: 921 1st Qu.:-200.0 1st Qu.: 4.005 1st Qu.: 711.0
## Median :1052 Median :-200.0 Median : 7.887 Median : 894.5
## Mean :1049 Mean :-159.1 Mean : 1.866 Mean : 894.5
## 3rd Qu.:1221 3rd Qu.:-200.0 3rd Qu.: 13.636 3rd Qu.:1104.8
## Max. :2040 Max. :1189.0 Max. : 63.741 Max. :2214.0
## NOx(GT) PT08.S3(NOx) NO2(GT) PT08.S4(NO2)
## Min. :-200.0 Min. :-200.0 Min. :-200.00 Min. :-200
## 1st Qu.: 50.0 1st Qu.: 637.0 1st Qu.: 53.00 1st Qu.:1185
## Median : 141.0 Median : 794.2 Median : 96.00 Median :1446
## Mean : 168.6 Mean : 794.9 Mean : 58.14 Mean :1391
## 3rd Qu.: 284.2 3rd Qu.: 960.2 3rd Qu.: 133.00 3rd Qu.:1662
## Max. :1479.0 Max. :2682.8 Max. : 339.70 Max. :2775
## PT08.S5(O3) T RH AH
## Min. :-200.0 Min. :-200.000 Min. :-200.00 Min. :-200.0000
## 1st Qu.: 699.8 1st Qu.: 10.950 1st Qu.: 34.05 1st Qu.: 0.6923
## Median : 942.0 Median : 17.200 Median : 48.55 Median : 0.9768
## Mean : 975.0 Mean : 9.777 Mean : 39.48 Mean : -6.8376
## 3rd Qu.:1255.2 3rd Qu.: 24.075 3rd Qu.: 61.88 3rd Qu.: 1.2962
## Max. :2522.8 Max. : 44.600 Max. : 88.72 Max. : 2.2310##Data Cleaning & Preprocessing Nilai -200 pada dataset merupakan kode untuk missing value, sehingga dikonversi menjadi NA. Selanjutnya hanya variabel numerik yang digunakan dalam analisis karena PCA dan FA berbasis matriks korelasi numerik.
#Ganti -200 menjadi NA
data[data == -200] <- NAAmbil variabel numerik
data_num <- data[, sapply(data, is.numeric)]Cek Persentase Missing Value
Persentase missing value dihitung untuk mengidentifikasi variabel dengan tingkat kehilangan data tinggi.
na_percent <- colSums(is.na(data_num)) / nrow(data_num) * 100
tabel_na <- data.frame(
Variabel = names(na_percent),
Persen_NA = round(na_percent, 2)
)
kable(tabel_na, caption = "Tabel Persentase Missing Value")| Variabel | Persen_NA | |
|---|---|---|
| CO(GT) | CO(GT) | 17.99 |
| PT08.S1(CO) | PT08.S1(CO) | 3.91 |
| NMHC(GT) | NMHC(GT) | 90.23 |
| C6H6(GT) | C6H6(GT) | 3.91 |
| PT08.S2(NMHC) | PT08.S2(NMHC) | 3.91 |
| NOx(GT) | NOx(GT) | 17.52 |
| PT08.S3(NOx) | PT08.S3(NOx) | 3.91 |
| NO2(GT) | NO2(GT) | 17.55 |
| PT08.S4(NO2) | PT08.S4(NO2) | 3.91 |
| PT08.S5(O3) | PT08.S5(O3) | 3.91 |
| T | T | 3.91 |
| RH | RH | 3.91 |
| AH | AH | 3.91 |
Drop NMHC(GT)
Variabel NMHC(GT) dihapus karena memiliki proporsi missing yang sangat besar dan berpotensi merusak struktur analisis.
data_num <- data_num[, colnames(data_num) != "NMHC(GT)"]Imputasi Median
Missing value yang tersisa diimputasi menggunakan median untuk menjaga kestabilan distribusi dan meminimalkan distorsi terhadap varians data.
data_num <- data_num %>%
mutate(across(everything(),
~ ifelse(is.na(.), median(., na.rm = TRUE), .)))Cek ulang NA
colSums(is.na(data_num))## 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
## RH AH
## 0 0##Karakteristik Data
Jumlah Observasi & Variabel
Menunjukkan banyaknya sampel dan jumlah variabel yang digunakan setelah proses pembersihan data.
nrow(data_num)## [1] 9357ncol(data_num)## [1] 12Statistik Deskriptif
untuk memahami sebaran dan skala masing-masing variabel sebelum dilakukan standardisasi.
desc_stats <- describe(data_num)
tabel_deskriptif <- data.frame(
Variabel = rownames(desc_stats),
Mean = round(desc_stats$mean, 3),
SD = round(desc_stats$sd, 3),
Min = round(desc_stats$min, 3),
Max = round(desc_stats$max, 3),
row.names = NULL
)
kable(tabel_deskriptif,
caption = "Tabel Statistik Deskriptif Variabel Kualitas Udara")| Variabel | Mean | SD | Min | Max |
|---|---|---|---|---|
| CO(GT) | 2.089 | 1.323 | 0.100 | 11.900 |
| PT08.S1(CO) | 1098.272 | 212.915 | 647.250 | 2039.750 |
| C6H6(GT) | 10.011 | 7.311 | 0.149 | 63.741 |
| PT08.S2(NMHC) | 937.855 | 261.623 | 383.250 | 2214.000 |
| NOx(GT) | 235.131 | 195.093 | 2.000 | 1479.000 |
| PT08.S3(NOx) | 834.203 | 251.808 | 322.000 | 2682.750 |
| NO2(GT) | 112.360 | 43.938 | 2.000 | 339.700 |
| PT08.S4(NO2) | 1456.402 | 339.368 | 551.000 | 2775.000 |
| PT08.S5(O3) | 1020.452 | 390.779 | 221.000 | 2522.750 |
| T | 18.294 | 8.659 | -1.900 | 44.600 |
| RH | 49.245 | 16.974 | 9.175 | 88.725 |
| AH | 1.024 | 0.396 | 0.185 | 2.231 |
#Visualisasi (Boxplot) Untuk melihat distribusi masing masing variabel, mengidentifikasi keberadaan outlier, dan membandingkan rentang nilai antar variabel, karena skala tiap variabel berbeda, standarisasi diperlukan sebelum PCA.
for(i in 1:ncol(data_num)){
boxplot(data_num[, i],
main = paste("Boxplot", colnames(data_num)[i]),
col = "lightblue",
ylab = colnames(data_num)[i])
}
Deteksi Outlier
Outlier dideteksi menggunakan metode IQR (Interquartile Range). Observasi yang berada di luar rentang; Q1−1.5×IQR dan Q3+1.5×IQR dianggap sebagai outlier. Outlier dapat memengaruhi struktur komponen utama karena PCA sensitif terhadap variansi ekstrem.
count_outliers <- function(x){
Q1 <- quantile(x, 0.25)
Q3 <- quantile(x, 0.75)
IQR <- Q3 - Q1
sum(x < (Q1 - 1.5*IQR) | x > (Q3 + 1.5*IQR))
}
sapply(data_num, count_outliers)## CO(GT) PT08.S1(CO) C6H6(GT) PT08.S2(NMHC) NOx(GT)
## 454 145 281 92 778
## PT08.S3(NOx) NO2(GT) PT08.S4(NO2) PT08.S5(O3) T
## 277 379 131 130 12
## RH AH
## 0 7Standardisasi Data
Data dinormalisasi menggunakan fungsi scale() agar memiliki: - Mean = 0 - Standar deviasi = 1
Standardisasi penting karena PCA berbasis kovarians/korelasi dan sangat dipengaruhi oleh skala variabel.
data_scaled <- scale(data_num)
head(data_scaled)## CO(GT) PT08.S1(CO) C6H6(GT) PT08.S2(NMHC) NOx(GT) PT08.S3(NOx)
## [1,] 0.38600794 1.2292598 0.2558843 0.411452879 -0.3543494 0.8818111
## [2,] -0.06749848 0.9110579 -0.0839442 0.064579235 -0.6772723 1.3484357
## [3,] 0.08367033 1.4265216 -0.1385655 0.005333599 -0.5337510 1.2144053
## [4,] 0.08367033 1.3020588 -0.1069730 0.039734291 -0.3235949 1.0237841
## [5,] -0.36983609 0.8171237 -0.4777149 -0.391229934 -0.5337510 1.4725380
## [6,] -0.67217370 0.4636965 -0.7207952 -0.717080932 -0.7490329 1.9947605
## NO2(GT) PT08.S4(NO2) PT08.S5(O3) T RH AH
## [1,] 0.01455834 0.69422651 0.6321924 -0.5420791 -0.02178476 -0.6734325
## [2,] -0.46338291 0.30158462 -0.1233488 -0.5767248 -0.09100664 -0.7549382
## [3,] 0.03731744 0.28906133 0.1370282 -0.7384045 0.27866730 -0.6924149
## [4,] 0.21939030 0.37525102 0.4677774 -0.8423415 0.63361313 -0.6002820
## [5,] 0.08283566 0.09900203 0.2291518 -0.8250187 0.60857550 -0.5950236
## [6,] -0.37234649 -0.18682358 -0.1822055 -0.8221315 0.58501060 -0.6051847##Uji Kelayakan PCA/FA Uji KMO mengukur kecukupan sampel untuk analisis faktor.
- KMO > 0.5 → Layak
- KMO > 0.7 → Baik
- KMO > 0.8 → Sangat baik
Uji Bartlett menguji apakah matriks korelasi berbeda signifikan dari matriks identitas.
- p-value < 0.05 → Variabel saling berkorelasi → Layak dilakukan PCA/FA
#KMO
KMO(data_scaled)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data_scaled)
## Overall MSA = 0.79
## MSA for each item =
## CO(GT) PT08.S1(CO) C6H6(GT) PT08.S2(NMHC) NOx(GT)
## 0.91 0.95 0.84 0.80 0.81
## PT08.S3(NOx) NO2(GT) PT08.S4(NO2) PT08.S5(O3) T
## 0.82 0.90 0.80 0.95 0.46
## RH AH
## 0.24 0.43#Bartlett Test
cortest.bartlett(cor(data_scaled), n = nrow(data_scaled))## $chisq
## [1] 164437.2
##
## $p.value
## [1] 0
##
## $df
## [1] 66Langkah perbaikan
Variabel dengan nilai MSA < 0.50 dihapus karena tidak cukup berkontribusi dalam struktur faktor.
Setelah penghapusan variabel (RH, T, AH), uji KMO dan Bartlett dilakukan kembali untuk memastikan kelayakan analisis meningkat.
data_reduced <- data_num[, !(colnames(data_num) %in% c("RH","T","AH"))]
data_scaled2 <- scale(data_reduced)Ulang KMO dan Bartlett Test (Setelah perbaikan)
#KMO
KMO(data_scaled2)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data_scaled2)
## Overall MSA = 0.87
## MSA for each item =
## CO(GT) PT08.S1(CO) C6H6(GT) PT08.S2(NMHC) NOx(GT)
## 0.91 0.94 0.82 0.80 0.85
## PT08.S3(NOx) NO2(GT) PT08.S4(NO2) PT08.S5(O3)
## 0.88 0.87 0.86 0.93#Bartlett
cortest.bartlett(cor(data_scaled2), n = nrow(data_scaled2))## $chisq
## [1] 120313.4
##
## $p.value
## [1] 0
##
## $df
## [1] 36##Principal Component Analysis (Unrotated) # PCA Full (untuk Eigenvalue) PCA pertama dilakukan tanpa rotasi untuk memperoleh eigenvalue dan menentukan jumlah komponen optimal.
pca_unrotated <- principal(data_scaled2,
nfactors = ncol(data_scaled2),
rotate = "none")eigenvalue
Jumlah komponen ditentukan berdasarkan: - Kriteria Kaiser (eigenvalue > 1) - Visualisasi scree plot (titik elbow)
Komponen sebelum titik elbow dianggap menjelaskan variansi terbesar secara signifikan.
eigen_values <- pca_unrotated$values
kable(eigen_values, caption = "Tabel Eigenvalue PCA")| x |
|---|
| 6.5746546 |
| 1.1999168 |
| 0.4306272 |
| 0.2454993 |
| 0.2018748 |
| 0.1352151 |
| 0.1167909 |
| 0.0838088 |
| 0.0116124 |
Scree Plot
plot(eigen_values, type="b",
main="Scree Plot",
xlab="Komponen",
ylab="Eigenvalue")
abline(h=1, col="red", lty=2)
PCA unrotated Final
PCA dijalankan kembali dengan 2 komponen utama sesuai hasil eigenvalue.
pca_2 <- principal(data_scaled2,
nfactors = 2,
rotate = "none")
print(pca_2)## Principal Components Analysis
## Call: principal(r = data_scaled2, nfactors = 2, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 h2 u2 com
## CO(GT) 0.88 0.15 0.80 0.201 1.1
## PT08.S1(CO) 0.93 -0.12 0.89 0.115 1.0
## C6H6(GT) 0.95 -0.20 0.94 0.065 1.1
## PT08.S2(NMHC) 0.96 -0.20 0.96 0.040 1.1
## NOx(GT) 0.76 0.54 0.87 0.128 1.8
## PT08.S3(NOx) -0.84 0.02 0.70 0.297 1.0
## NO2(GT) 0.70 0.60 0.85 0.153 2.0
## PT08.S4(NO2) 0.69 -0.65 0.90 0.101 2.0
## PT08.S5(O3) 0.94 0.01 0.88 0.125 1.0
##
## PC1 PC2
## SS loadings 6.57 1.20
## Proportion Var 0.73 0.13
## Cumulative Var 0.73 0.86
## Proportion Explained 0.85 0.15
## Cumulative Proportion 0.85 1.00
##
## Mean item complexity = 1.3
## Test of the hypothesis that 2 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.03
## with the empirical chi square 524.05 with prob < 5e-99
##
## Fit based upon off diagonal values = 1PCA dengan Rotasi Varimax
(Ganti nfactors sesuai eigenvalue > 1) Rotasi Varimax dilakukan untuk: - Memperjelas struktur loading - Menghasilkan interpretasi yang lebih mudah - Mengurangi kompleksitas item
pca_rot <- principal(data_scaled2,
nfactors = 2,
rotate = "varimax")
print(pca_rot)## Principal Components Analysis
## Call: principal(r = data_scaled2, nfactors = 2, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
## RC1 RC2 h2 u2 com
## CO(GT) 0.59 0.68 0.80 0.201 2.0
## PT08.S1(CO) 0.80 0.50 0.89 0.115 1.7
## C6H6(GT) 0.86 0.44 0.94 0.065 1.5
## PT08.S2(NMHC) 0.87 0.45 0.96 0.040 1.5
## NOx(GT) 0.24 0.90 0.87 0.128 1.1
## PT08.S3(NOx) -0.66 -0.52 0.70 0.297 1.9
## NO2(GT) 0.16 0.91 0.85 0.153 1.1
## PT08.S4(NO2) 0.95 -0.06 0.90 0.101 1.0
## PT08.S5(O3) 0.72 0.60 0.88 0.125 1.9
##
## RC1 RC2
## SS loadings 4.40 3.37
## Proportion Var 0.49 0.37
## Cumulative Var 0.49 0.86
## Proportion Explained 0.57 0.43
## Cumulative Proportion 0.57 1.00
##
## Mean item complexity = 1.5
## Test of the hypothesis that 2 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.03
## with the empirical chi square 524.05 with prob < 5e-99
##
## Fit based upon off diagonal values = 1Komunalitas
Komunalitas menunjukkan proporsi variansi masing-masing variabel yang dapat dijelaskan oleh komponen yang dipilih.
Nilai mendekati 1 → variabel sangat terwakili oleh komponen. Nilai rendah → variabel kurang cocok dalam model.
kable(pca_rot$communality,
caption = "Tabel Komunalitas PCA (Rotasi Varimax)")| x | |
|---|---|
| CO(GT) | 0.7992868 |
| PT08.S1(CO) | 0.8850289 |
| C6H6(GT) | 0.9352207 |
| PT08.S2(NMHC) | 0.9597948 |
| NOx(GT) | 0.8716150 |
| PT08.S3(NOx) | 0.7025209 |
| NO2(GT) | 0.8466049 |
| PT08.S4(NO2) | 0.8991830 |
| PT08.S5(O3) | 0.8753164 |
##Factor Analysis (FA) FA bertujuan mengidentifikasi faktor laten yang mendasari korelasi antar variabel. # Menentukan Jumlah Faktor (Parallel Analysis) Parallel analysis digunakan untuk menentukan jumlah faktor optimal dengan membandingkan eigenvalue aktual dengan eigenvalue acak.
Faktor dipilih jika eigenvalue aktual > eigenvalue acak.
fa.parallel(data_scaled2, fa = "fa", fm = "ml")
## Parallel analysis suggests that the number of factors = 3 and the number of components = NAJalankan FA
Sesuaikan jumlah faktor hasil parallel
fa_result <- fa(data_scaled2,
nfactors = 2,
rotate = "varimax",
fm = "ml")
print(fa_result$loadings, cutoff = 0.4)##
## Loadings:
## ML1 ML2
## CO(GT) 0.539 0.698
## PT08.S1(CO) 0.738 0.526
## C6H6(GT) 0.865 0.469
## PT08.S2(NMHC) 0.883 0.464
## NOx(GT) 0.917
## PT08.S3(NOx) -0.631 -0.500
## NO2(GT) 0.795
## PT08.S4(NO2) 0.871
## PT08.S5(O3) 0.676 0.614
##
## ML1 ML2
## SS loadings 4.064 3.299
## Proportion Var 0.452 0.367
## Cumulative Var 0.452 0.818Komunalitas FA
Menunjukkan seberapa besar variansi tiap variabel dijelaskan oleh faktor laten.
Nilai tinggi menunjukkan model faktor mampu merepresentasikan variabel dengan baik.
kable(fa_result$communality,
caption = "Tabel Komunalitas FA")| x | |
|---|---|
| CO(GT) | 0.7771272 |
| PT08.S1(CO) | 0.8221636 |
| C6H6(GT) | 0.9683431 |
| PT08.S2(NMHC) | 0.9942634 |
| NOx(GT) | 0.8840978 |
| PT08.S3(NOx) | 0.6476375 |
| NO2(GT) | 0.6772710 |
| PT08.S4(NO2) | 0.7583875 |
| PT08.S5(O3) | 0.8335953 |
Visualisasi Loading
Heatmap digunakan untuk: - Memvisualisasikan kekuatan hubungan variabel terhadap faktor - Mengidentifikasi pola klaster variabel - Mempermudah interpretasi struktur faktor
par(mar = c(1, 1, 3, 1))
corrplot(as.matrix(fa_result$loadings),
is.corr = FALSE,
tl.cex = 0.7)
title("Heatmap Factor Loadings", line = 1)