Importing Data

#Importing library
library("readxl")
library("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, union
library("MVN")
## Warning: package 'MVN' was built under R version 4.3.3
library("psych")
## Warning: package 'psych' was built under R version 4.3.3
#Importing Data
df <- read_excel("Data APG_Pulau Jawa_Kelompok 3.xlsx")

#Statistik Deskriptif
xbar <- colMeans(df[,2:5])
s <- cov(df[,2:5])

Uji Independensi

cortest.bartlett(cor(df[,2:5]), n = 119)
## $chisq
## [1] 49.04731
## 
## $p.value
## [1] 7.295155e-09
## 
## $df
## [1] 6
qchisq(p = 1-0.05, df = (4+2)*(4-1)/2)
## [1] 16.91898

Hipotesis

\[ H_0:\Sigma=\sigma^2I~(Antar~variabel~independen)\\ H_1:\Sigma\neq\sigma^2I~(Antar~variabel~tidak~independen) \]

\[ \alpha=5\% \]

\[ statistik~uji=-2[1-(2p^2+p+2)/6pn]~ln~\Lambda \]

\[ Tolak~H_0~jika~statistik~uji>\chi^2_{(p+2)(p-1)/2;\alpha} \]

\[ statistik~uji = 49.047\\ titik~kritis=\chi^2_{(4+2)(4-1)/2;0.05}=16.92\\ Keputusan:Tolak~H_0 \]

Pengecekkan Outlier

Chisquare Plot

#Pembuatan Chisquare Plot
d2 <- mahalanobis(df[,2:5], xbar, s)
sorted_d2 <- sort(d2)
Q <- rep(0, 119)
for (i in 1:119) {
  Q[i] <- qchisq(p = (i-0.5)/119, df = 2)
}

#Evaluasi Grafis Menggunakan Chisquare Plot
plot(x = Q, 
     y = sorted_d2,
     xlim = c(0,6.5),
     ylim = c(0,6.5))
abline(a=0,b=1, col="red")

cor(Q, sorted_d2)
## [1] 0.744623

Interpetasi:

Berdasarkan grafik yang diberikan, terdapat beberapa titik yang melenceng dari garis y = x. Hal ini mengindikasikan bahwa data yang dimiliki terindikasi memiliki observasi yang bersifat outlier. Namun, untuk lebih memberikan gambaran yang jelas terkait dengan keberadaan outlier, maka akan dilakukan proses deteksi outlier dengan metode mahalanobis distance.

Deteksi Outlier (Mahalanobis)

Berikut adalah alur kalkulasi yang digunakan untuk deteksi outlier pada metode Mahalanobis:

  1. Hitung Mahalanobis distance dari data menggunakan rumus berikut:

    \[ d_i^2 = (x_i - \vec{\mu})^T\Sigma^{-1}(x_i - \vec{\mu}) \]

  2. Bandingkan nilai distance terhadap titik kritis dari \(\chi^2_{1-\alpha,df~=~p}\) dengan p adalah banyak variabel atau dimensi dari data. Jika \(d_i^2>\chi_{1-\alpha;df~=~p}^2\) maka data tersebut merupakan outlier.

#Penghitungan titik kritis
cval <- qchisq(p = 1-0.05, df = 4)

#Penentuan Outlier
outlier <- ifelse(d2 > cval, FALSE, TRUE)

#Penentuan Baris Outlier
df_cleaned <- df[outlier,]
outliers <- df[!outlier,]

outliers
## # A tibble: 8 × 5
##   kab         jmh_penduduk_ribu persen_angkatan_kerja   rls rata_pendapatan_juta
##   <chr>                   <dbl>                 <dbl> <dbl>                <dbl>
## 1 Kota Tange…             1951.                  93.1 10.9                  4.04
## 2 Cimahi                   591.                 895.  11.4                  2.39
## 3 Jakarta Se…             2236.                  94.6 11.8                  4.45
## 4 Jakarta Ut…             1779.                  93.0 10.8                  4.75
## 5 Banjarnega…             1047.                 937.   6.86                 1.24
## 6 Rembang                66017                   97.4  7.72                 1.78
## 7 Kabupaten …             1007.                 968.   7.47                 1.67
## 8 Kabupaten …             1177.                 968.   6.29                 1.46

Interpretasi:

Berdasarkan data tersebut, terdapat 7 outlier yakni sesuai dengan table yang tertera. Berikutnya, akan dilakukan proses PCA.

Analisis Komponen Utama

Misal:

\[ X1: Jumlah~penduduk~dalam~ribuan\\ X2: Persentase~penduduk~angkatan~kerja\\ X3: Rata-rata~lama~sekolah\\ X4: Rata-rata~pendapatan~dalam ~Juta~Rupiah \]

Maka, berikut ini adalah proses perhitungan PCA-nya:

PCA tanpa eliminasi outliers

#PCA Covariance
pca_cov <- prcomp(df[,2:5], center = TRUE ,scale. = FALSE)
pca_cov$rotation
##                                 PC1           PC2          PC3           PC4
## jmh_penduduk_ribu      9.999998e-01 -0.0005563989 1.618225e-05  9.318516e-06
## persen_angkatan_kerja -5.564129e-04 -0.9999994618 8.756172e-04  6.696915e-06
## rls                   -1.758681e-05  0.0008468154 9.650920e-01  2.619098e-01
## rata_pendapatan_juta   4.878944e-06  0.0002228669 2.619097e-01 -9.650923e-01
pca_cov$sdev
## [1] 5994.9543370  155.6803358    1.6674737    0.5333705
#PCA Correlation
pca_cor <- prcomp(df[,2:5], scale. = TRUE)
pca_cor$rotation
##                               PC1        PC2         PC3         PC4
## jmh_penduduk_ribu     -0.02052505  0.8590123 -0.49455638 -0.13073125
## persen_angkatan_kerja -0.15542670 -0.4996759 -0.85128532 -0.03846834
## rls                    0.70114618 -0.1089452 -0.03225876 -0.70390649
## rata_pendapatan_juta   0.69556833  0.0235129 -0.17229817  0.69709768
pca_cor$sdev
## [1] 1.2611541 1.0108753 0.9869034 0.6431511

PCA Matriks Kovarians

Berikut adalah komponen utama jika menggunakan matriks kovarians:

\[ Y_1=0.999X_1-0.000556X_2-0.00001759X3+0.00000488X4\\ Y_2=-0.000556X_1-0.99999X_2+0.000847X_3+0.000223X_4\\ Y_3=0.0000162X_1+0.000876X_2+0.965X_3+0.262X_4\\ Y_4=0.00000932X_1+0.000006697X_2+0.262X_3-0.965X_4 \]

Berikut adalah Varians yang dijelaskan oleh masing-masing komponen utama:

Komponen Utama Sdev Varians
Y1 5994.954 35939473.462
Y2 155.68 24236.2624
Y3 1.667 2.779
Y4 0.533 0.284

PCA Matriks Korelasi

Berikut adalah komponen utama jika menggunakan matriks kovarians:

\[ Y_1=-0.0205X_1-0.155X_2+0.701X_3+0.696X_4\\ Y_2=0.859X_1-0.500X_2-0.109X_3+0.023X_4\\ Y_3=-0.495X_1-0.851X_2-0.0322X_3-0.172X_4\\ Y_4=-0.131X_1-0.038X_2-0.704X_3+0.697X_4 \]

Berikut adalah Varians yang dijelaskan oleh masing-masing komponen utama:

Komponen Utama Cumulative Varians
Y1 0.25
Y2 0.50
Y3 0.75
Y4 1.00

PCA dengan Eliminasi Outliers

#PCA Covariance without outliers
pca_cov_cleaned <- prcomp(df_cleaned[,2:5], center = TRUE, scale. = FALSE)
pca_cov_cleaned$rotation
##                                 PC1           PC2           PC3           PC4
## jmh_penduduk_ribu     -0.9999983483 -0.0017898053  0.0001314032  0.0002874992
## persen_angkatan_kerja -0.0017907714  0.9999891261 -0.0040021666 -0.0015885478
## rls                    0.0001824472  0.0042502064  0.9766791376  0.2146619774
## rata_pendapatan_juta  -0.0002513492 -0.0006928395  0.2146667825 -0.9766870683
pca_cov_cleaned$sdev
## [1] 890.5636978  20.7642135   1.6010251   0.3747661
#PCA Correlation
pca_cor_cleaned <- prcomp(df_cleaned[,2:5], center = TRUE, scale. = TRUE)
round(pca_cor_cleaned$rotation,4)
##                           PC1     PC2     PC3     PC4
## jmh_penduduk_ribu     -0.3840 -0.7655 -0.2298  0.4624
## persen_angkatan_kerja -0.0938 -0.2867  0.9500 -0.0805
## rls                   -0.5630  0.5760  0.1664  0.5688
## rata_pendapatan_juta  -0.7257 -0.0048 -0.1304 -0.6755
pca_cor_cleaned$sdev
## [1] 1.2794036 1.0529763 0.9945613 0.5149904

PCA Matriks Kovarians

Berikut adalah komponen utama jika menggunakan matriks kovarians:

\[ Y_1=-0.9999X_1-0.00179X_2+0.00018X_3-0.000251X_4\\ Y_2=-0.000179X_1+0.9999X_2+0.00425X_3-0.000693X_4\\ Y_3=0.000131X_1-0.004X_2+0.977X_3+0.215X_4\\ Y_4=0.0003X_1-0.0016X_2+0.215X_3-0.977X_4 \]

Berikut adalah varians dari komponen utama yang diperoleh:

Komponen Utama Sdev Varians
Y1 890.56 793097.1136
Y2 20.76 430.9776
Y3 1.601 2.563
Y4 0.3747 0.1404

PCA Matriks Korelasi

Berikut adalah komponen utama jika menggunakan matriks korelasi:

\[ Y_1=-0.384X_1-0.0938X_2-0.563X_3-0.726X_4\\ Y_2=-0.765X_1-0.287X_2+0.576X_3-0.0048X_4\\ Y_3=-0.230X_1+0.950X_2+0.166X_3-0.130X_4\\ Y_4=0.462X_1-0.0805X_2+0.569X_3-0.675X_4 \]

Berikut adalah varians dari komponen utama yang diperoleh:

Komponen Utama Cumulative Varians
Y1 0.25
Y2 0.50
Y3 0.75
Y4 1.00

Scree Plot

#dengan outlier
##PCA Matriks kovarians
plot(pca_cov, type = "lines", main = "Scree Plot PCA Covariance dengan Outlier")

##PCA Matriks korelasi
plot(pca_cor, type = "lines", main = "Scree Plot PCA Corellation dengan Outlier")

#Tanpa Outlier
##PCA Matriks kovarians
plot(pca_cov_cleaned, type = "lines", main = "Scree Plot PCA Covariance tanpa Outlier")

##PCA Matriks korelasi
plot(pca_cor_cleaned, type = "lines", main = "Scree Plot PCA Corellation tanpa  Outlier")